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A  high  resolution  frequency-wavenumber  spectral  estimation 
technique  is  employed  to  analyze  spatio-temporal  signals  from  three- 
dimensional  sensor  arrays,  to  estimate  the  spectrum  for  finite 
temporal  frequency  bands  using  FFT  techniques,  and  to  estimate  the 
spectrum  for  discrete  temporal  frequencies. 

Direction  and  true  vector  velocity  of  a  plane  wavefront  propa- 
gating across  a  three-dimensional  arrayare  given  by  the  frequency- 
wavenumber  spectrum  for  discrete  temporal  frequencies.   For  finite 
temporal  frequency  bands  the  wavefront  direction  and  bounds  on  the 
vector  velocity  are  given. 

A  parametric  analysis  is  performed  in  which  the  properties  of 
the  frequency-wavenumber  spectral  estimate  are  determined  for  known 
signals  and  parametric  values.  Array  aperture  and  the  ratio  of 
incoherent  noise  power  to  total  signal  plus  noise  power  are  found  to 
determine  the  resolution  in  the  wavenumber  domain.  Other  parameters, 
including  number  of  temporal  samples  per  sensor,  temporal  window  function, 
and  number  of  sensors,  are  found  to  affect  statistical  reliability  and 
computation  cost. 


The  frequency-wavenumber  spectral  estimation  technique  is 
applied  to  the  analysis  of  visual-evoked  response  (VER)  data  and  to 
the  analysis  of  penicillin-induced  epileptiform  electrical  activity 
recorded  from  the  rat  neocortex.  Results  indicate  that  application 
of  the  frequency-wavenumber  spectral  analysis  method  to  electro- 
physiological data  is  possible  and  that  it  provides  information  for 
increased  understanding  of  brain  function. 

Application  of  the  expanded  frequency-wavenumber  spectral 
estimation  method  to  spatio-temporal  data  from  the  areas  of  seismology, 
oceanography,  speech,  meteorology,  radar  and  sonar  would  require  in 
most  cases  only  a  change  in  parametric  inputs. 
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CHAPTER  I 
INTRODUCTION 

Over  the  past  few  years  Biomedical  Engineering  has  evolved 
from  isolated  applications  of  electronic  instruments  as  tools  for 
medicine  to  sophisticated  and  ever-expanding  applications  of  advanced 
analytical  techniques  to  biological  systems.  Other  areas  of  applica- 
tion, which  include  seismology,  oceanography,  speech,  meteorology 
and  sonar,  have  need  for  analytical  techniques  similar  to  those 
applicable  to  the  analysis  of  biological  data.   These  various 
disciplines  have  certain  common  characteristics  which  makes  it 
advantageous  to  develop  analytical  techniques  that  have  general 
applicability.  One  such  characteristic  is  that  each  of  the  areas  is 
concerned  with  measuring  and  interpreting  spatio-temporal  signals. 
The  purpose  of  this  research  was  to  investigate  a  technique  for 
analyzing  spatio-temporal  signals  and  develop  useful  extensions  of 
the  technique  to  provide  additional  information  about  the  signals. 
Parametric  bounds  and  example  computations  were  chosen  to  emphasize 
application  of  the  technique  to  the  study  of  electrophysiological 
data.  However,  application  to  other  areas  is  implied  and  in  most 
cases  would  require  only  a  change  in  parametric  inputs. 

Chapter  11  gives  a  brief  review  of  array-, processing  techniques 
which  have  been  applied  to  other  fields  such  as  radar,  seismology, 
acoustics,  and  image  processing.  Some  of  the  properties  of  these 


techniques  are  discussed  and  compared  with  the  expected  properties  of 
the  frequency-wavenumber  technique.  Digital  filtering  array  processors 
are  indicated  as  having  general  applicability  as  wavenumber  limiting 
pre-f liters. 

In  Chapter  III  the  theoretical  development  is  given  for  the 
frequency-wavenumber  spectral  estimate.   It  is  shown  how  the  frequency- 
wavenumber  spectrum  provides  information  for  determining  direction  and 
velocity  for  multiple  plane  wavefronts.  An  extension  of  the  frequency- 
wavenumber  technique  to  three  dimensions  is  developed  and  shown  to  be 
valid  theoretically.   Inclusion  of  frequency  bandpass  signals  as  well 
as  discrete  frequency  signals  has  been  accomplished  using  fast  Fourier 
transform  (FFT)  techniques  and  a  spectral  Integration  scheme.  Details 
of  the  bandpass  feature  are  given  in  Chapter  III. 

Many  parameters  affect  the  resolution,  cost  and  reliability 
of  the  frequency-wavenumber  spectral  estimate.  Additive  incoherent 
noise  was  found  to  have  a  major  effect  on  resolution;  however,  the 
array  aperture  and  the  temporal  window  were  found  to  interact  with 
the  incoherent  noise  level  and  thus  affected  the  achievable  resolution. 
Parametric  variation,  statistical  reliability  and  computation  require- 
ments are  discussed  in  Chapter  IV. 

For  spherically  sjnmnetric  wavefronts  the  frequency-wavenumber 
method  gives  a  scalar  output  which  is  insufficient  for  determining 
wavefront  direction.  Chapter  IV  discusses  the  case  for  spherical 
wavefronts  as  well  as  the  need  for  spatial  pre-f ilterlng. 

In  Chapter  V  examples  are  shown  for  simulation  runs  which 
verify  the  theoretical  development  and  parametric  analyses . 


The  frequency-wavenumber  technique  was  applied  to  electro- 
physiological data  of  two  types:  1.  hmnan  visual-evoked  response 
data  recorded  from  the  scalp  over  the  occipital  region  of  the  brain, 
and  2.  penicillin-induced  focal  epileptic  discharges  recorded  from 
the  rat  neocortex.  Details  and  results  of  these  experiments  are 
outlined  in  Chapter  VI. 

Finally  in  Chapter  VII,  a  discussion  of  the  frequency-wavenumber 
spectral  estimation  technique  is  given  in  relation  to  the  goals  accom- 
plished by  this  research  effort.  Significance  of  the  experimental 
results  is  discussed  along  with  recommended  areas  for  further  investiga- 
tions. 


CHAPTER  II 
ARRAY  PROCESSING 

Included  under  the  category  of  array  processing  are  several 
specific  techniques  which  have  been  investigated  by  various  researchers. 
The  characteristics  of  some  of   these  techniques  are  discussed  in  rela- 
tion to  the  array-processing  technique  termed  frequency-wavenumber 
spectral  estimation. 

In  general,   array- processing  techniques  are  used  in  applications 
where  their  advantages  outweigh  the  disadvantage  of  increased  processing 
complexity.      The  advantages  offered  by  array- processing  techniques 
include  improved  signal-to-noise  ratio  for  facilitating  the  detection  of 
weak  signals,  enhancement  or  emphasis  of   specific  features  or  parameters 
and  the  determination  of  spatial  characteristics  such  as  wavefront  direc- 
tion.    Major  differences  between  the  various  array- processing  techniques 
usually  result  from:      1.  noise   assumptions    (stationarity ,  directionality); 
2.   signal  assumptions    (known,   unknown,  uniform  or  distorted  across  the 
array);   and  3.   method  employed    (spatial   filtering,    stochastic  approxima- 
tion, steepest  descent,   conjugate  gradients,  noise  estimation,  etc.). 
Various  methods  of  array-processing  are  discussed  below. 

A.     Adaptive  Processors  for  Signal-to-Noise  Ratio  Improvement 

Although  many  techniques  have  been  developed  for  improving 
signal-to-noise  ratio  using  array  data,   most  of  the   techniques  are 


similar  in  approach  and  make  the  assumption  that  wavefront  direction  is 
known  "a  priori." 

Lacoss  [1]  has  presented  an  adaptive  linear  processor  with 
variable  coefficients  in  which  the  coefficients  are  adjusted  by  a  rule 
similar  to  that  used  by  the  projection  gradient  method  of  a  quadratic 
form  subject  to  a  linear  constraint.  The  purpose  of  this  processor  is 
to  improve  signal  detection  and  assumes  that  the  wavefront  direction  is 
known  with  'sufficient'  accuracy.  Thus  the  signals  from  the  sensors 
can  be  combined  after  appropriate  time  delay  adjustments  with  gains 
adjusted  to  provide  optimum  detection  of  an  unknown  signal.  An  optimum 
processor  is  one  which  converges  to  a  minimum  variance,  unbiased  estimate 
of  the  signal.  Lacoss  derived  several  iterative  methods  for  reducing 
sensitivity  of  the  processor  to  data  anomalies  such  as  noise  bursts  and 
for  reducing  required  computation  time  and  storage  requirements.  Other 
adaptive  processing  techniques  for  improving  slgnal-to-noise  ratio  which 
are  similar  to  those  of  Lacoss  have  been  derived  in  the  literature. 
Shor  [2]  has  proposed  a  gradient  adaptive  method  for  narrowband  hydrophone 
arrays  which  requires  knowledge  of  the  signals'  autocorrelation  function. 
Adams  [3]  has  done  work  with  adaptive  on-line  array  processors  for 
receiving  deep-space  probe  signals. 

Kobayashi  [4]  has  proposed  two  Iterative  methods  for  processing 
array  data:  the  method  of  steepest  descent  and  the  method  of  conjugate 
gradients  with  projection.  These  methods  are  very  similar  to  those  of 
Lacoss  and  rely  on  essentially  the  same  mathematical  development.  No 
intermediate  computation  of  covariance  matrix  functions  is  required 
for  the  methods  used  by  Kobayashi.  As  do  most  methods,  however,  this 


method  assumes  the  only  difference  between  the  signals  at  various  sensors 
is  a  time  delay  which  can  be  compensated. 

Burg  [5]  used  a  multi-dimensional  Wiener  filtering  approach  to 
array  processing.   In  this  approach  the  signal  and  noise  are  represented 
as  stationary  random  processes  with  known  cross-correlation  functions. 
The  output  of  the  processor  is  a  minimum  mean-squared-error  estimate  of 
the  signal.  Again,  knowledge  of  time  delays  between  sensors  is  required. 

'  Another  approach  to  array  processing  is  that  proposed  by  Claerbout 
[6],  Tlie  cross-correlation  matrix  is  computed  by  observing  the  sensor 
outputs  in  a  "fitting  interval."  This  information  is  then  used  to  provide 
a  minimum  mean-squared-error  estimate  of  the  noise  for  a  short,  future, 
projected  interval.  The  predicted  noise  is  then  subtracted  from  the 
signal.   This  method  reduces  the  noise  level  but  introduces  distortion 
into  the  signal. 

Capon  et^al^.  [7]  have  discussed  the  advantages  and  disadvantages 
of  time-domain  versus  frequency-domain  array  processors.   The  signal-to- 
noise  improvement  is  best  with  time-domain  methods;  however,  time-domain 
methods  are  more  sensitive  to  non-stationary  noise  and  signal  anomalies 
than  the  frequency-domain  methods. 

B.   Feature  Extraction  by  Digital  Filtering 

The  fields  of  feature  extraction  and  digital  filtering  are  by 
far  too  extensive  to  be  given  a  comprehensive  survey  in  this  document; 
however,  feature  extraction  from  two-dimensional  fields  by  digital 
filtering  techniques  is  a  form  of  array  processing  and  is  included  for 
that  reason.   In  addition  to  feature  extraction,  digital  filters  may 


be  used  in  array  processors  to  limit  wavenuniber  content  and  reduce 
aliasing  effects   for  a  given  spatial  sampling  rate. 

Wavenumber  limited  spatial  filters  are  discussed  further  in 
Chapter  IV. 

C.     Wavefront  Detection 

To  determine  a  wavefront  propagation  direction  it  is  necessary 
to  use  methods  which  give  vector  wavenumbers  as  outputs.      Standard 
beamforming  techniques  use  adjustable  phase  differences  between 
elements   of  a  sensor  array  to  determine  wavefront  directions.      These 
methods,  however,   provide  no  information  regarding  wave  velocities  and 
are  restricted  in  resolution  to   the  natural  beam  pattern  of  the  array. 
A  method  which  overcomes  these  deficiencies   and  computes  the  high 
resolution  vector  velocity  and  direction  of  propagating  wavefronts 
has  been  proposed  by  Capon   [8]    for  processing  two-dimensional  data 
from  the  large-aperture  seismic  array   (LASA)   in  Montana.     A  comparison 
of  this  high  resolution  estimate  of  the  frequency-wavenumber  spectrum 
with  conventional  estimates  showed  marked  improvement   for  the  high 
resolution  method.      The  theoretical  development  of  the   frequency- 
wavenumber  spectrum  is  given  in  Chapter  III.      Extension  of  Capon's 
method  to  include  data  from  three-dimensional  arrays   is  shown  to  be 
valid  theoretically.     Methods   for  including  finite  frequency  bands 
are  discussed  and  one  method  is  developed  in  detail. 


CHAPTER  III 
FREQUENCY-WAVENUMBER  SPECTRUM—THEORY 

This  chapter  gives  a  detailed  development  of  the  theoretical 
background  material  supporting  the  validity  of  the  frequency-wavenumber 
technique  as  a  descriptor  of  array  sensed  signals.  First  the  properties 
of  plane  waves  and  their  mathematical  representations  are  discussed 
and  it  is  shown  how  the  frequency-wavenumber  spectrum  provides  informa- 
tion describing  plane  waves. 

The  concept  of  estimation  as  related  to  the  frequency-wavenumber 
spectrum  is  introduced  and  related  to  theoretical  concepts.  Window 
functions  are  discussed  in  relation  to  their  effect  on  accuracy  of  the 
estimate.  A  high  resolution  adaptive  window  is  introduced  which  provides 
a  major  improvement  in  resolution  over  that  obtainable  with  conventional 
fixed  window  functions. 

The  use  of  a  two-dimensional  array  for  measuring  three-dimensional 
signals  is  inadequate  if  it  is  desired  to  know  true  vector  velocities  of 
incoming  wavefronts.   Determination  of  wavenumbers  in  two  dimensions 
identifies  azimuth  angle  but  not  elevation  angle.  This  means  that  for  a 
given  azimuth  angle,  any  number  of  signals  with  different  elevation  angles 
may  be  present.   With  the  two-dimensional  analysis  these  various  signals 
would  be  indistinguishable  and  only  their  phase  velocities  could  be  deter- 
mined.  A  three-dimensional  analysis  identifies  wavenumbers  in  three 
dimensions  and  uniquely  determines  wavefront  directions.   This  also  makes 
it  possible  to  determine  true  vector  velocity. 


Theory  and  methods  are  discussed  for  extending  the  frequency- 
wavenumber  technique  to  accept  data  from  a  three-dimensional  array  of 
sensors  and  to  accept  finite  bandwidth  data  in  addition  to  discrete 
frequencies.  Capon's  technique  is  valid  for  only  single  discrete 
frequencies.  This  could  be  a  major  disadvantage  where  the  exact  spectral 
content  of  the  signal  is  unknown.  Also  in  many  applications  it  may  be 
desirable  to  include  a  band  of  frequencies  even  if  spectral  content  is 
known.  Rather  than  compute  spectral  estimates  individually  for  all 
frequencies  of  interest  a  method  was  developed  for  accepting  multiple 
frequencies  as  well  as  finite  continuous  frequency  bands  and  efficiently 
computing  the  frequency-wavenumber  spectrum. 

A.  Generalized  Plane  Waves 

A  general  plane  wave  in  two  dimensions  can  be  represented 
mathematically  by   [9], [10] 

S(t,x,y)  =  exp  [2Tri(ft  -  k  -r)]  (1) 

where:       S(t,x,y)  =  the  spatio-temporal  signal  amplitude 

f  =  temporal  frequency  (Hz) 
k  =  vector  wavenumber  (cycles /m) 
r  =  spatial  coordinate  vector  (m) 

For  non-isotropic  media,  the  vector  wavenumber  (k)  is  spatially 
dependent  and  is  a  function  of  x  and  y.   For  an  attenuating  medium  k 
has  both  real  and  imaginary  components.  Thus  the  general  vector  wave- 
number  is  given  by 

1  =  Re[k  (x,y)]  +  Im[k  (x,y)]  (2) 

S  S  a 
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where:       Re[']  denotes  the  "real  part  of" 

Im[*]  denotes  the  "imaginary  part  of" 

If  the  signal  is  represented  in  an  isotropic,  non-attenuating 

medium  then  k  is  a  real,  vector  constant  and  is  denoted  by  k.   The 
g 

spatio-temporal  signal  representation  for  a  two-dimensional,  isotropic, 
non-attenuating  medium  then  becomes 

S(t,x,y)  =  exp  [2TTi(ft  -  k-r)]  (3) 

Figure  1  shows  the  above  relationships  and  how  the  vector  wave- 
number  specifies  the  propagation  direction  for  the  plane  wavefront. 
For  simplicity  and  clarity  in  demonstrating  the  frequency-wavenumber 
technique  an  isotropic,  non-attenuating  medium  will  be  assumed  so  that 
for  the  following  analysis  equation  (3)  is  representative  of  a  plane 
wave.   Some  of  the  properties  of  a  plane  wave  or  any  general  spatio- 
temporal  signal  which  is  to  be  measured  by  an  array  of  sensors  are 
often  more  easily  studied  by  measuring  the  cross  power  spectrum  between 
any  two  sensors.   The  cross  power  spectral  density  is  a  K  x  K  matrix 
for  an  array  of  K  sensors  and  is  related  to  the  coherency  matrix  which 
provides  a  measure  of  the  linear  dependence  among  the  sensors.   The 
cross  power  spectral  density  matrix  may  be  obtained  in  either  of  two 
ways:   the  correlation  method  or  the  direct  method.   Details  of  these 
two  methods  are  given  in  Appendix  A. 

Table  1  provides  a  comparison  of  pertinent  criteria  for  the 
correlation  method  versus  the  direct  method  for  estimating  the  cross 
power  spectral  density.   Because  of  reduced  computation  time  and  less 
susceptibility  to  signal  variations  the  direct  method  was  chosen  for 
the  estimate.   Statistical  reliability  is  improved  by  segmenting  the 
data  as  described  by  Welch  [11].   In  some  cases  however,  increased 
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>x 


ikr  =  3-^  +  7' 

k       =  vector  wavenumber   (gives  wavefront  direction) 


GENERAL : 

S(t,x,y) 

k 
g 


exp   [27ri(ft  -  k  'r)] 
Re[kg(x,y)]   +  Im[kg(x,y)] 
X    +  y 


ISOTROPIC,    NON-ATTENUATING  MEDIA: 

S(t,x,y)   =  exp    [2iri(ft  -  k*r)] 

=  exp  [2iri(ft  -  gx  -  Yy)] 
WHERE:  k  =  Re[k  (x,y)]  =  constant 

3  =  component  of  k  in  x  direction 
Y  =  component  of  k  in  y  direction 


Figure  1.   GENERAL  TWO-DIMENSIONAL  PLANE  WAVE 
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statistical  reliability  is  achieved  at  the  cost  of  reduced  resolving 
capability. 


Table  1:   CROSS  POWER  SPECTRAL  DENSITY  ESTIMATE 


CORRELATION  METHOD 

DIRECT  METHOD 

DEFINITION 

J^     k=-L/2  J^ 

L  =  NO.  OF  DATA  POINTS 
p,„(k)  IS  COVARIANCE 
MATRIX 

F,  (A)  IS  FOURIER  TRANSFORM 
OF  DATA  IN  mth  SEGMENT  OF 
jth  SENSOR. 

M  SEGMENTS  OF  N  DATA  POINTS 
EACH  FOR  EACH  SENSOR 

COMPUTATIONAL 
TIME  REQUIRED 

MORE 

LESS 

STATISTICAL 
RELIABILITY 

GOOD  FOR  LARGE  L 

GOOD  FOR  LARGE  M  &  N 

SUSCEPTIBILITY 

TO  SIGNAL  VARIA- 
TION 

MORE 

LESS 

B.  Frequency-Wavenumber  Representation  of  Spatial  Fields 

In  Appendix  A  it  is  shown  how  the  cross  power  spectral  density 
function  for  K  sensors  is  a  K  x  K  matrix.   The  cross  power  spectral 
density  matrix  is  a  temporal  frequency  transform  of  the  cross  covariance 
matrix  for  an  array  and  its  measured  signals.   The  cross  power  spectral 
density  matrix  contains  information  concerning  the  frequency  content 
of  the  signal  power.   In  a  similar  manner  it  is  possible  to  obtain 
information  about  the  wavenumber  content  of  the  signal  by  performing  a 
spatial  frequency  transform  on  the  cross  power  spectral  density  function. 
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The  result  is  a  frequency-wavenumber  spectrum  and  is  defined  as 


00 

pa,k)  =    fa,s)  e^^^^'^  dx  (4) 

—00 

Vector 


where  f(X,k)  is  the  cross  power  spectral  density  function  and  the 
integral  is  a  vector  Fourier  transform  with  respect  to  the  spatial 
vector,  X.   It  follows  that  the  cross  power  spectral  density  is  then 
the  Fourier  inverse  of  the  frequency-wavenumber  spectrum. 


fa,x)  = 


F(X,k)  e"^^^^*''  dk  (5) 


—00 

Vector 


This  development  is  based  on  a  treatment  by  Yaglom  [12]  and  details 
of  the  development  are  given  in  Appendix  B. 

It  is  also  shown  in  Appendix  B  that  if  the  signal  waveform  is 
a  discrete  frequency  plane  wave  then  the  frequency-wavenumber  spectrum 
is  a  multi-dimensional  delta  function  about  the  given  frequency  and 
vector  wavenumber ,  i.e., 

P(X,k)  =  6(X  -  X^)  6(k  -  k^)  (6) 

Thus  it  is  shown  that  wavefront  detection  is  possible  through  the  use 
of  the  frequency-wavenumber  spectrum. 

C.   Estimate  of  the  Frequency -Wavenumber  Spectrum 

An  exact  spectral  representation  for  spatio-temporal  signals 
would  require  infinite  summations  over  time  and  space.  Since  this  is 
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not  practical  a  spectral  representation  using  a  finite  number  of  sensors 
(spatial  representation)  and  a  finite  number  of  time  samples  (temporal 
representation)  is  computed  and  distinguished  from  the  exact  representa- 
tion by  being  called  an  estimate  of  the  exact  spectrum.   Criteria  for 
measuring  goodness  of  an  estimate  have  been  developed  and  are  discussed 
in  Chapter  IV,  Section  E  on  statistical  reliability.  As  was  shown 
previously  there  are  two  major  methods  for  estimating  power  spectral 
densities,  the  correlation  method  and  the  direct  method.   Because  of 
its  desirable  statistical  characteristics  and  computational  efficiency 
the  direct  segment  method  or  block  averaging  technique  was  used  for 
estimating  the  cross  power  spectral  density.   The  resulting  matrix  is 
denoted  bylJ.5(X)]  to  distinguish  it  as  an  estimate  of  the  exact  spectrum 
f  .(X).   This  estimate  of  the  cross  power  spectral  density  is  then  used 
to  obtain  an  estimate  of  the  frequency-wavenumber  spectrum. 

To  illustrate  the  technique  assume  that  discrete  data  are  avail- 
able from  K  sensors.  Data  from  each  of  the  sensors  aire  divided  into  M 
non-overlapping  segments  of  N  data  points  each  to  give  L  =  MN  data  points 
per  sensor.  The  cross  power  spectral  density  matrix  is  computed  using 
the  direct  segment  method.  The  Fourier  transform  of  the  data  in  the  mth 
segment  of  the  jth  sensor  is 


F,  (X)  =  -T7T  I     ^n  S,„(nT)  e^"^  (7) 

j  =  1, ,K 

m  =  1, ,M 
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where:  S.  (nT)  =  the  data  in  the  mth  segment  from  the  jth  sensor 
X  =  ZirfT  =  normalized  frequency 

a  =  the  weighting  coefficients  describing  the 
temporal  window  used  in  estimating  the 
frequency  spectrum. 

For  simplicity  assume  a  =  1;  n  =  1, ,N.   In  practice,  however,  side- 
lobe  leakage  can  be  reduced  and  the  estimate  improved  in  some  cases 
by  using  windows  different  from  the  rectangular  window  described  here. 

From  (7)  an  estimate  for  the  cross  power  spectral  density  is 
a  K  X  K  matrix  whose  elements  are  given  by 


•'  m=l  •' 


In  order  .to  remove  the  effect  of  differences  in  gain  and  other  causes 
of  improper  sensor  equalization,  a  normalization  is  performed  on  f^^C^) 

The  cross  power  spectral  density  matrix  is  closely  related  to 
coherency,  the  major  difference  being  that  phase  information  is  retained 
in  the  normalized  cross  power  spectrum  and  lost  in  the  coherency  function. 
The  relationship  is  given  by 

2 

where  Yjo^^)  ^^  ^^^   coherency  function.   The  coherency  function  varies 
J*' 
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from  zero  to  one  and  is  a  measure  of  linear  dependence  between  sensor 
elements . 

Using  the  normalized  cross  power  spectral  density  estimate  and 
the  coordinates  of  the  sensors,  we  can  obtain  the  conventional  estimate 
for  the  frequency-wavenumber  spectrum: 


,        K       K  ik2iT'(x.-  Xj) 

Pa,k)  =K     1       I     w.w*  f.<,(A)  e  J  (11) 


'a,k)  =  ^     I       I     w.w*  f  .(A) 


where  the  w.  are  the  weights  describing  the  shape  of  the  spatial  window 
used  to  estimate  the  wavenumber  spectrum.  Again  for  simplicity  we  assume 

w  =1;  j  =  1, ,K  so  that  w,  represents  the  aperture  of  the  sensor 

array.  As  was  indicated  for  the  temporal  window,  the  spatial  window 
can  be  designed  to  provide  an  improved  estimate  of  the  spectrum.  One 
such  window,  proposed  by  Capon  [8],  is  a  function  of  the  particular 
wavenumber  being  considered.   It  is  different  for  each  wavenumber  k 
and  passes  undistorted  any  monochromatic  plane  wave  traveling  with  a 
velocity  corresponding  to  the  wavenumber  k  .   It  also  suppresses  in 
an  optimum  least  squares  sense  the  power  of  waves  traveling  at  veloc- 
ities  corresponding  to  wavenumbers  other  than  k  . 

Figure  2  shows  a  flow  chart  describing  the  estimation  procedure 
for  obtaining  the  frequency-wavenumber  spectrum.   It  provides  a  summary 
of  the  procedure  discussed  up  to  this  point  and  the  resultant  output 
is  called  the  conventional  frequency-wavenumber  spectrum  of  a  homogeneous 
random  field. 
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Figure  2.      CONVENTIONAL  FREQUENCY-WAVENUMBER  SPECTRAL  ESTIMATE 
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D.  High  Resolution  Adaptive  Window 

A  finite  data  record  of  discrete  values  may  be  considered  as 
the  product  of  an  Infinite  representation  of  the  data  with  a  rectangular 
sampling  window  as  shown  In  Figure  3a.  From  Fourier  transform  theory 
we  know  that  multiplication  in  the  time  domain  Implies  convolution  In 
the  frequency  domain.   The  result  of  the  convolution  Is  shown  In 
Figure  3b  for  the  time  functions  in  3a.  An  ideal  situation  would  be 
to  have  a  window  function  whose  transform  is  an  impulse.   This  of 
course  Implies  a  window  function  of  infinite  width  which  implies  an 
infinite  data  record  length.   Thus  it  is  desirable  to  find  a  window 
function  which  more  closely  approximates  the  ideal  but  with  a  finite 
data  record.  More  detail  on  window  functions  and  their  characteristics 
which  have  been  used  in  this  research  is  given  under  Section  C  of 
Chapter  IV  on  the  parametric  effects  of  temporal  window  design.  The 
choice  of  which  window  to  use  depends  on  the  particular  application 
and  criteria  selected  to  determine  a  desirable  output.  There  is  no 
window  function  which  is  generally  accepted  as  optimum. 

Capon  Introduces  a  high  resolution  estimate  for  the  frequency- 
wavenumber  spectrum  which  utilizes  a  window  function  that  is  dependent 
on  the  frequency  and  wavenumber  being  considered.  The  high  resolution 
frequency-wavenumber  spectral  estimate  using  this  window  is  given  by 

K   K  12 Trie*  Cx   ~  x  ) 

P'a,k)=  I       I     A*a,k)  A  (A,k)  f.j(X)  e      ^   ^    (12) 
j=l  Jl=l  J       ^       ^^ 

where:       A,  (A,k)  Is  the  high  resolution  window  function. 
Comparison  of  equations  (11)  and  (12)  Illustrates  the  fact  that  the 
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high  resolution  estimate  is  the  direct  result  of  choosing  an  adaptive 

window  function.   The  window  function  is  given  by 

K 

Aj  a  ,k)  =  -Y^ (13) 

j=l  £=1  J* 

where:        [qjj^a.k)]  =  [  [exp  {ZT^lt- (k^   -  Xj^)} ] •  [fj^(X)]  ]~^    (14) 


that  is,  [q.(,]  is  the  inverse  of  the  matrix  resulting  from  the  product 
of  the  exponential  matrix  [exp  (27rik»  (x  -  Xp))]  with  the  cross  power 


t' 


spectral  density  matrix  [fjn] 


Further,  if  [q.n]  is  defined  as  the  inverse  of  [f^nl 


[q.,]  =  [fj,]-'  (15) 


then  it  can  be  shown  that  an  equivalent  expression  for  the  high  resolu- 
tion frequency-wavenumber  spectral  estimate  is  given  by 


P'(X,k)  = 


K        K 

I       I     q.„(A)  exp    [2Trik-(x     -  x  ): 

J=l  £=1     ^^  ^         ^ 


(16) 


Thus  P'  can  be  computed  by  simply  taking  the  Hermitian  inverse  of  the 
cross  power  spectral  density  matrix  and  then  computing  the  spatial 
transform.  Equivalency  of  the  two  expressions,  (12)  and  (16),  is 
difficult  to  show  in  general;  however,  computations  with  low  order 
matrices  have  verified  their  validity. 

The  estimate  P'  is  the  power  output  of  a  maximum  likelihood 
array  processor  whose  characteristics  vary  with  wavenumber  in  an 
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Figure  4.   HIGH  RESOLUTION  FREQUENCY-WAVENUMBER  SPECTRAL  ESTIMATE 
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optimum  least  squares  sense.   That  is,  P'  provides  a  minimum  variance, 
unbiased  estimate  of  the  power.  For  each  wavenumber,  k  ,  it  passes 
undistorted  any  monochromatic  plane  wave  with  a  velocity  corresponding 
to  the  wavenumber  k  and  suppresses  in  an  optimum  least  squares  sense 
all  waves  traveling  at  velocities  other  than  k  . 

•Figure  4  gives  a  flow  chart  outlining  the  procedure  for  obtaining 
the  high  resolution  frequency-wavenumber  spectral  estimate.   The  major 
difference  between  it  and  the  conventional  technique  is  the  inclusion 
of  the  matrix  inversion  step.   Since  [f.g]   is  in  general  complex, 
additional  care  must  be  exercised  in  implementing  the  technique. 

E.  Three-Dimensionsl  Frequency-Wavenumber  Spectrum 

The  theoretical  development  of  the  frequency-wavenumber  spectrum 
is  not  limited  to  a  specific  number  of  dimensions.   Therefore  the  equations 
developed  for  the  two-dimensional  frequency-wavenumber  spectrum  are 
directly  applicable  to  detection  of  three-dimensional  wavefronts  if  the 
vectors  k  and  x  are  extended  to  include  components  in  three  dimensions. 
The  detection  of  wavefronts  in  three  dimensions  would  be  particularly 
useful  for  the  measurement  of  seismological  events  via  earth  body  waves 
and  the  measurement  of  EEG  signals  which  can  certainly  be  considered  to 
possess  three-dimensional  characteristics.  Although  the  theoretical 
development  is  straightforward,  difficulties  arise  in  the  implementation 
of  a  three-dimensional  technique.  A  pictorial  representation  of  three- 
dimensional  wavefront  characteristics  is  given  below.  Then  some  of  the 
implementation  difficulties  are  discussed  further. 

Figure  5  shows  how  a  plane  wave  might  be  represented  in  three 

dimensions.  Determination  of  the  three  wavenumber  components  k  ,  k  , 

X   y 
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Figure  5.  THREE-DIMENSIONAL  PLANE  WAVE 
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and  k.  specifies  not  only  the  azimuth  and  phase  velocity  of  the  signal 
but  provides  the  added  advantage  of  specifying  the  elevation  angle 
and  group  velocity  of  the  signal.  Determination  of  the  wavenuraber  com- 
ponents is  accomplished  through  the  estimation  of  the  three-dimensional 
frequency-wavenuniber  spectrum. 

The  major  problems  involved  In  the  extension  of  the  frequency- 
wavenumber  spectrum  to  three  dimensions  are  complexity  of  implementation, 
increased  cost,  how  to  display  the  data  in  a  meaningful  way  and  how  to 
design  sensor  arrays  to  best  sense  data  in  specific  applications  such 
as  biological  preparations.  An  order  of  magnitude  increase  in  computa- 
tion time  and  thus  cost  results  from  computation  of  the  spatial  transform 
in  three  dimensions  over  that  required  for  computation  in  two  dimensions. 
Fast  Fourier  transform  (FFT)  techniques  are  not  readily  applicable  to 
the  spatial  transform  operation.   Since  most  FFT  algorithms  require  that 
the  number  of  data  points  be  a  power  of  two  and  since  all  FFT  techniques 
are  based  on  data  points  at  equal  sampling  intervals,  the  use  of  these 
techniques  places  stringent  requirements  on  the  sensor  arrays  to  be  used. 
Equal  spacing  of  electrodes  for  the  measurement  of  biological  signals 
may  in  some  cases  be  impossible  or  even  undesirable.   If  this  is  so  then 
a  more  general  technique  for  determining  the  frequency-wavenumber  spectrum 
is  required,  which  precludes  use  of  the  fast  Fourier  transform. 

Extension  of  the  frequency-wavenumber  technique  to  three  dimensions 
was  accomplished  and  test  cases  were  run  using  both  known  input  signals 
and  data  recorded  from  the  rat  neocortex  using  a  three-dimensional  electrode 
array.  Details  of  the  results  of  the  three-dimensional  frequency-wavenumber 
analysis  are  discussed  in  Chapter  V  for  the  known  input  signal. 
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Meaningful  display  methods  for  the  three-dimensional  spectrum 
depend  to  some  extent  on  the  wavefront  complexity.  A  method  which 
provides  two-dimensional  spectra  for  sequential  values  of  the  third 
dimension  was  used  for  all  three-dimensional  data  runs  because  of  its 
ease  of  comparison  with  the  two-dimensional  data  runs.  This  method 
worked  quite  well  for  simple,  known  input  signals;  however,  for  more 
complicated  signals  it  proved  to  be  less  clear  for  demonstrating 
wavefront  directions. 

The  particular  array  geometry  to  be  used  for  the  three- 
dimensional  spectral  estimate  is  strongly  dependent  on  the  specific 
application.  However,  to  demonstrate  many  of  the  principles  and 
parametric  effects,  a  simple  cubic  array  geometry  with  known  input 
signals  was  used.  For  practical  application  of  the  technique  to 
remotely  sensed  spatio-temporal  signals  such  as  the  EEG,  certain 
assumptions  about  the  characteristics  of  the  signal  were  made. 
These  assumptions  are  discussed  in  Chapter  VI  along  with  experimental 
results. 

F.  Finite  Bandpass  Signals 

Recall  from  the  equations  defining  the  conventional  and  high 
resolution  frequency-wavenumber  spectral  estimates  that  both  are  func- 
tions of  A  the  normalized  frequency.  This  says  that  for  each  frequency 
contained  in  the  signal  a  different  wavenumber  spectrum  can  be  computed. 
Clearly  this  procedure  is  not  practical  for  the  case  where  exact 
frequency  content  is  not  known  or  where  a  large  range  of  frequency 
components  is  known  to  compose  the  signal.  The  question  then  arises 
as  to  how  the  technique  may  be  modified  to  accept  a  narrow  finite  band 
of  frequencies. 
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Consider  the  case  of  a  composite  signal  given  by 


s  =  5!  s  (t,x,y,z)  (17) 


where:  s  (t,x,y,z)  =  exp  [2Trl(f.t  -k  x-k  y-k  z)]        (18) 
^  ^     ""j     ^j     ^j 

In  order  for  the  technique  to  identify  the  components  of  the  signal, 
it  must  compute  the  cross  power  spectral  density  estimate  which  includes 
all  the  frequencies ,  f . .  For  a  small  number  of  frequencies  this  may 
be  accomplished  by  computing  each  component  and  then  adding  to  get  the 
total.  This  particular  technique  was  tested  for  the  case  of  two 
frequency  components  and  appeared  to  introduce  very  little  distortion 
into  the  frequency-wavenumber  spectral  estimate.  Up  to  five  discrete 
frequencies  have  been  included  in  the  spectral  estimation  technique 
with  only  a  minor  increase  in  computation  time.  For  large  numbers 
of  discrete  components  or  where  interest  in  finite  bandpass  signals 
exists  an  integration  technique  on  the  output  from  an  FFT  algorithm 
can  be  used.  The  inclusion  of  multiple  frequencies  introduces  some 
distortion  into  the  estimation  of  the  frequency-wavenumber  spectrum, 
as  does  the  finite  bandpass  technique.  The  extent  of  this  distortion 
is  determined  by  a  complex  interaction  of  parameters  including  incoherent 
noise  ratio,  data  length,  and  temporal  window  function. 

A  method  for  including  finite  bandpass  data  was  developed  which 
uses  the  average  integrated  spectral  output  over  the  selected  spectral 
band.  The  IBM  subroutine,  HARM,  was  used  in  the  fast  Fourier  transform 
spectral  computations .   If  the  output  of  HARM  is  as  illustrated  in 
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Figure  6  where  the  spectrum  shovm  Is  for  the  jth  sensor,  mth  data 
segment ,  then  an  estimate  for  the  bandpass  spectral  content  is 

IHIGH 
F.^a)  =    I       A  -Af  (19) 

-^       i=ILOW 

where:       A.  =  spectral  value  at  ith  frequency 

FLOW 


ILOW  = 


Af 


IHIGH  =  ^IP-  +  1 

FLOW  =  low  frequency  cuton 
FHIGH  =  high  frequency  cutoff 

Af  =  frequency  resolution  =  r— r ;— -; r- 

record  length  (sec) 

The  Fjm(^)  are  then  used  in  the  direct  segment  method  for  computing 
cross  power  spectral  density. 

It  is  important  to  note  that  the  inclusion  of  multiple  frequencies 
either  as  discrete  components  or  finite  bands  introduces  ambiguity  into 
determination  of  phase  velocities  or  group  velocities.  However,  if 
wavefront  directions  are  of  prime  importance  then  this  method  of  accepting 
multiple  frequencies  will  provide  more  information  more  efficiently. 
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CHAPTER  IV 
FREQUENCY-WAVENUMBER  SPECTRUM—PARAMETRIC  ANALYSIS 

In  general  any  spectral  analysis  technique  should  produce  an 
estimate  that  is  statistically  reliable  and  which  has  high  resolution. 
It  is  of  interest  therefore  to  identify  parameters  affecting  the 
estimate  and  to  quantify  their  effects  on  resolution,  statistical 
reliability,  and  other  measures  of  goodness  of  the  estimate. 

The  spatial  and  temporal  sample  intervals.  Ax  and  AT,  are 
determined  primarily  by  the  expected  bandwidths  in  their  respective 
transform  domains.   These  sample  intervals  are  chosen  to  reduce 
aliasing  in  the  transform  domain  of  the  periodically  repeated 
spectrum.   Temporal  resolution  is  shovm  to  be  inversely  proportional 
to  record  length  and  spatial  resolution  is  inversely  proportional  to 
array  aperture.  Therefore  for  fixed  sample  intervals.  Ax  and  AT, 
resolution  is  improved  by  increasing  the  number  of  sample  values. 
The  effect  of  an  increased  number  of  sample  values  on  cost  is  partic- 
ularly important  when  discrete  signals  in  time  and  three  spatial 
coordinates  are  considered  such  as  those  for  the  three-dimensional 
frequency-wavenumber  spectral  estimate. 

Another  parameter  which  results  from  choosing  the  direct  segment 
method  for  estimating  the  cross  power  spectral  density  matrix  is  the 
number  of  segments  M  into  which  each  data  record  is  divided.  The 
value  of  M  has  an  effect  on  the  statistical  reliability  of  the  frequency- 
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wavenumber  spectral  estimate  and  should  be  large ,  which  implies  more 
data  points  and  greater  cost  so  a  trade-off  is  indicated. 

Temporal  window  design  can  be  utilized  to  reduce  sidelobe 
leakage  and  in  some  cases  improve  resolution.  Various  window  func- 
tions and  their  characteristics  are  discussed. 

Included  under  the  parametric  analysis  are  discussions  of 
spherical  wavefronts  in  relation  to  the  objective  of  determining 
wavefront  direction  and  spatial  pre-filters  in  relation  to  expected 
maximum  wavenumbers. 

A.   Comparison  of  High  Resolution  and  Conventional  Estimates 

It  was  stated  that  Capon's  high  resolution  spatial  window 
function  improves  resolution  over  that  obtainable  with  a  conventional 
window  function.   It  is  therefore  of  interest  to  give  a  quantitative 
comparison  of  the  conventional  and  high  resolution  estimates  in  rela- 
tion to  achievable  resolution.  The  following  analysis  is  derived 
largely  from  Capon's  paper  [8], 

Consider  a  single  plane  wave  (wavenumber  =  k  )  propagating 
across  an  array  of  sensors  plus  a  component  of  noise  in  each  sensor 
which  is  incoherent  between  any  pair  of  sensors.  Assume  that  M,  the 
number  of  data  segments  per  sensor,  and  N,  the  number  of  data  points 
per  segment,  are  large  so  that  an  unbiased,  consistent  estimate  of 
the  frequency-wavenumber  spectrum  is  obtained.  Then  the  cross  power 
spectral  density  matrix  elements  are  estimated  as 


f.j^(A^)  =  6.^(R)  exp  [-i2TTk^«(Xj  -  Xj^)]   j  ,£  =  1, ,K   (20) 
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where:       A  =  2'irf  T 
o     o 


and 


f  =  temporal  frequency  of  incident  wave  (Hz) 

C  =  vector  wavenumber  of  incident  wave  (cycles/m) 

X.  =  vector  position  coordinates  of  jth  sensor 


^*  (21) 


=  1-R  ;     j  ^  il 


„ incoherent  noise  power  ^^o-v 

total  power  of  sensor  output 


Using  the  above  expression  for  the  input  signal  and  the  equations  for 
conventional  and  high  resolution  estimates  of  the  frequency-wavenumber 
spectrum  it  can  be  shown  that 


P(X^,S)  =  (l-R)|B(Ak^)|^  +  R/K  (23) 

and 


(24) 
(25) 

B(k)  =  i  ^  exp  [2Trik'x  ]  (26) 

j=l  ^ 


and  K  equals  the  number  of  sensors.   For  k  =  k 


P(X^,k^)  =  P'(Ajj,k^)  =  1-R  +  R/K  (27) 


°              "^  1-R  +  2  (R/K)   -  P(X^) 

where : 

Aic    =  k  -  k 
0                    o 

and 
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In  the  neighborhood  of  k   ,    consider  those  values  of  k  such 

o 

that  the  conventional  estimate  gives  a  value  close  to  Its  peak  value 

of  P(X  ,k  ). 
o'  o 

Since  R  Is  small   (0  ^  R  i  1)   and  K  Is  usually  large   (about  16) 
the  neighborhood  for  the  conventional  estimate  can  be  quite  large. 
For  these  values  of  k,    the  high  resolution  estimate  P'(X   ,k)   becomes 


^a^,k;       ^  ^_^  _^  2(R/K)    -   (1-R)  ^^^> 


or 


P'(A^,k)   =  I  (1-R+  R/K)  (29) 


which  is  already  3  db  down  from  its  peak  value.  Thus  P'  has  a  much 
steeper  falloff  and  gives  a  higher  resolution  estimate  for  the 
frequency-wavenumber  spectrum. 

B.  Array  Geometry 

The  data  to  be  represented  by  their  frequency-wavenumber  spectrum 
consist  of  finite  groups  of  spatio-temporal  samples.  That  is,  spatial 
sampling  is  a  result  of  a  finite  number  of  discrete  sensors  located 
within  a  finite  aperture.  Temporal  sampling  is  the  result  of  discrete 
sampling  of  a  finite  time  record  from  each  sensor.   The  relationship 
between  these  parameters  and  the  frequency-wavenumber  spectrum  is  given 
by  an  application  of  the  sampling  theorem,  modified  to  include  multi- 
dimensional signals.   For  example  the  transform  pair  represented  in 
Figure  7  demonstrates  properties  which  hold  for  any  bandlimited, 
truncated,  discrete  signal  which  includes  temporal  and  spatial  signals. 
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REQUIRE:   AT  <  1/2B 
RESOLUTION  IN  3  PLANE  =  1/T 

Figure  7.   SAMPLING  THEOREM  RELATIONSHIPS 

In  general,  the  spectral  resolution  varies  inversely  with  the  record 
length  and  the  required  sample  interval  varies  inversely  with  the 
bandwidth  of  the  signal. 

Sampling  theorem  concepts  can  be  applied  to  determine  the 
effects  of  array  geometry  on  the  resolution  of  the  frequency-wavenumber 
spectrum.  Application  of  these  concepts  is,  however,  complicated  by 
the  additional  property  of  dimensionality  for  spatial  fields.  For 
example,  equal  sample  intervals  for  a  spatial  field  has  meaning  only 
if  referenced  to  some  direction  in  the  field.   For  a  three-dimensional 
field,  equal  sampling  intervals  can  be  achieved  in  directions  represented 
by  three  perpendicular  coordinate  axes  by  constructing  a  cubic  sensor 
array.  However,  given  no  physical  limitations  on  the  array  size  or 
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shape,  Petersen  and  Middleton  [13]  have  shown  that  the  optimum  array 
(i.e. ,  minimum  number  of  sample  points  per  unit  hypervolume)  is  not 
in  general  rectangular.  The  optimum  array  is  represented  by  the 
centers  of  hyperspheres  in  their  closest  packed  arrangement.  For 
two  dimensions  the  optimum  array  is  rhombic  and  for  three  dimensions 
the  tightest  packing  of  spheres  is  achieved  with  a  hexagonal-close- 
packed  configuration.  These  arrangements  not  only  provide  maximum 
sampling  rates  for  a  given  number  of  electrodes  but  also  provide 
equal  sample  intervals  between  any  two  adjacent  sample  points.  For 
an  equal  number  of  elements  and  a  given  spatial  sampling  rate,  arrays 
other  than  these  have  a  reduced  array  aperture  and  resolution  is  not 
as  good. 

For  any  particular  array  design  the  sampling  theorem  relation- 
ship still  holds  that  resolution  in  a  given  wavenumber  component  is 
Inversely  proportional  to  the  aperture  of  the  array  in  that  dimension. 
In  most  applications  the  array  design  is  detemilned  by  physical 
limitations  of  the  recording  area.  This  is  particularly  true  in  the 
case  of  human  EEG  recordings  in  which  the  array  geometry  is  determined 
by  scalp  contour.   The  rationale  for  array  geometries  used  in  the 
experimental  data  examples  is  given  in  Section  G  and  in  Chapter  VI. 

C.   Temporal  Window  Functions 

It  was  shown  in  Chapter  III  that  the  spectrum  of  a  truncated 
signal  is  given  by  the  convolution  of  the  signal  transform  with  the 
transform  of  the  window  function.  An  ideal  temporal  window  function 
would  require  an  Infinite  length  of  time ;  therefore ,  it  is  of  interest 
to  investigate  various  finite  temporal  window  fimctlons  relative  to 
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their  sldelobe  levels,  bandwidths,  sharpness  of  cutoff,  and  statistical 

reliability.  The  simplest  window  function  is  a  rectangle  which  has  a 

sin  X 
transform  of  the  form  .   In  general  no  window  function  is  optimum. 

Rather,  the  best  window  design  for  a  particular  application  depends  on 
the  requirements  for  resolution,  statistical  reliability  and  sldelobe 
leakage. 

Several  window  functions  and  their  characteristics  are  shown 
in  Figure  8.  For  the  frequency-wavenumber  analyses  two  window  functions 
were  used.  The  rectangular  window  function  was  chosen  for  its  simplicity 
and  as  a  basis  for  comparison.  The  cosine  tapered  window  function  was 
also  used.   Because  of  its  reduced  sldelobe  level  the  cosine  window  is 
expected  to  provide  a  more  reliable  estimate  than  the  rectangular  window. 
For  the  frequency-wavenumber  spectral  estimate,  the  particular  temporal 
window  used  had  little  effect  on  spectral  resolution  in  comparison  with 
other  parameters  such  as  additive  incoherent  noise  level  and  array 
aperture.  Effects  of  the  various  parameters  on  resolution  are  demon- 
strated in  Chapter  V. 

D.  Additive  Incoherent  Noise 

Slncel  the  high  resolution  estimate  requires  a  matrix  Inversion 
of  the  cross  power  spectral  density  matrix  [f],  it  is  of  Importance 
to  determine  under  what  conditions  this  inverse  exists. 

From  matrix  theory  It  is  known  that  an  Hermltlan  matrix  has 
an  inverse  if  the  matrix  is  positive  definite  [15].   Capon  et  al.  [7] 
have  shown  that  the  cross  power  spectral  density  matrix  [f]  when 
computed  by  the  block  averaging  method  is  non-negative  definite. 
However,  this  is  not  sufficient  to  guarantee  that  an  Inverse  exists. 
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In  fact,  if  the  number  of  data  blocks  M  is  less  than  the  number  of 
sensors  K  then  the  matrix  is  singular.  This  comes  from  the  fact 
that  [f ]  is  the  sum  of  M  matrices  of  rank  unity. 


m=l 


where:  [q   ]   is  a  1  x  K  row  matrix  whose  jth  element  is  given  by 

_     r     „  inA 

*^nij        J^^  ''j.n+N(m-l)    ^  (31) 

Since  the  rank  of  [f ]  cannot  exceed  the  sum  of  the  ranks  of  its 
constituent  parts,  then  [f]  is  a  K  order  matrix  with  rank  M.  Clearly 
if  M  <  K  the  matrix  is  singular.  Thus  a  necessary  but  still  not 
sufficient  condition  for  [B]  to  have  an  inverse  is  that  M  ^  K. 

To  insure  that  an  inverse  exists  the  matrix  [f ]  is  modified 
to  give  the  matrix  f '  which  can  be  shown  to  be  positive  definite. 

[f]  =  (1  -  R)[f]  +  R[I]  (32) 

where:    R  =  ratio  of  incoherent  noise  power  to  total  power  at  a 

sensor 
and:      [I]  is  the  identity  matrix 

Thus  [f]  is  a  K  X  K  matrix  with  elements  given  by 


f  •  =  (1-R)  f   +  R  ;    m  =  n 

mm  ^33j 

=  (1-R)  f^^     ;    m  ^  n 
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Now  [f]  is  positive  definite  if  for  any  function  a.,  the  quadratic 
form  Q  given  by 


K   K 

j=l  Z=l     ^  ^   J^ 


(34) 


is  positive.  Using  the  definition  for  [f]  and  [f],  we  have 


Q  = 


1=R  V 


K   N 


inX 


K 


-  TT  <  X  ^  TT   (35) 

The  first  term  above  is  always  non-negative  for  it  represents 
a  positive  constant  (1-R)  multiplied  by  the  quadratic  form  for  [J] 
which  has  previously  been  demonstrated  to  be  non-negative  definite. 
Thus  the  only  way  for  Q  to  be  zero  is  if 


K 


(36) 


But  this  implies  that  a,  =  0;  j  =  1, ,K  which  is  a  contradiction. 

Therefore  [f]  is  positive  definite  and  its  inverse  exists. 

The  quantitative  dependence  of  the  high  resolution  frequency- 
wavenumber  estimate  on  the  magnitude  of  additive  incoherent  noise  is 
derived  from  expressions  given  in  Section  A  and  is  given  by 


P'(^o'^>  =f 


1-R  +  (R/K) 


1-R  +  (R/K) 

K   2Trl(It-E  )'x 


(1-R) 


^A^ 


(37) 
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Recall  that  the  peak  value  for  P'(A  ,k=k  )  is 


P'a^.l^o)  =  1-R+  (R/K) 


(38) 


so  that 


P'(A„,K) 

P'(A  ,Ic  ) 
^  o*  o 


(R/K)  +  (1-R) 


R/K 


1  - 


1   K   27ii(E-E  )'x.  1 
lie 
j=l 


(39) 


Effects  on  resolution  of  the  additive  incoherent  noise  are 
demonstrated  in  Chapter  V  by  artificially  injecting  noise  into  a 
known  signal.    (  See  also  Appendix  F.) 

E.   Statistical  Reliability 

Two  desirable  properties  of  any  estimator  are  that  it  be  un- 
biased and  consistent.  For  the  random  process  P  it  is  then  required 
of  the  estimate  P.,  that 


E[Pn^  =  ^ 


(unbiased) 


(40) 


lim  E[(P  -  f  )'']  =  0   (consistent) 


That  is,  the  mean  value  of  the  estimate  should  equal  the  mean  value  of 
the  random  process  and  as  the  number  of  obseirvations  goes  to  infinity, 
the  mean  square  error  between  the  estimate  and  the  random  process  should 
go  to  zero. 

Capon  and  Goodman  [16]  show  that  the  conventional  and  high  resolu- 
tion frequency-wavenumber  spectral  estimates  are  asymptotically  unbiased. 
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that  Is,  In  the  limit  as  N  (number  of  data  points  per  segment)  gets 
large  the  bias  approaches  zero.  The  expected  value  for  the  conven- 
tional estimate  Is  given  by 

E[P(X,S)]=|  I  P(x,S)  |B(k-k^)|2.|Wj^(x-A)|2|^dk        (41) 

-TT  -00 

Vector 

K       127rk»x 
where:  B(k)    =  ^    T     e  ^  (42) 

^  j=l 

Is  the  beamformlng  response  pattern  of  the  array  and. 


K(x)|'=| 


sin  (N/2)x 


(43) 


Is  the  Bartlett  window  resulting  from  the  rectangular  time  window  chosen 
for  the  temporal  transform. 

Thus  P  will  be  an  unbiased  estimator  for  P  If  In  the  limit  the 
frequency-wavenumber  window  approaches  a  delta  function, 

TT   oo 

I  I  |Wjj(x-X)|2.|B(k-ko^|^||cik^  6(x-A^)  6(k-k^)        (44) 

-TT  -00 

Vector 

Following  Goodman  [17]  and  Capon  and  Goodman  [16],  It  can  be 
shown  using  the  Complex  Wlshart  distribution  that  the  high  resolution 
estimate  is  unbiased  if  it  is  multiplied  by  M/(M-K+1).  That  is  the 
expected  value  of  P'  is 
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ir  °° 


E[P»a,k)]  =^^|  JPCX.IJ)  |B'a.k,k^)|2.|w^(x-X)|2|^dE 


-TT  -00 

Vector  (45) 


K 

where:     B'(X,k,k^)  =  I     A  (A,k^)  exp  [27ri(k-k^)'x  ]       (46) 


is  the  high  resolution  beamforming  pattern  of  the  array  and  is  a  func- 
tion of  k  .  Since  in  many  cases  all  spectral  values  are  normalized, 
the  question  of  bias  is  not  so  important.  Also  for  M  »  K  the  high 
resolution  estimate  gives  the  same  expected  value  as  the  conventional 
estimate. 

Assuming  that  the  noise  received  by  the  array  is  a  multi- 
dimensional Gaussian  process,  it  can  be  shown  for  large  N  that  the 
variance  of  the  conventional  estimate  is  given  by 

VAR  [P(X,k^)]  =  I  {E[P(A,k^)]}2  ;    |kj  =  0 

(47) 

=^{E[P(X,k  )]}2  ;    |k  I  /o 


Since  the  variance  goes  to  zero  as  M  goes  to  infinity  the  estimate  is 
consistent.  Again  using  the  Complex  Wishart  distribution  it  can  be 
shown  that  the  high  resolution  estimate  is  a  multiple  of  a  chi-square 
variable  with  2(M-K+1)  degrees  of  freedom  and  variance  given  by 

VAR  [P'(X,k^)]  =  jj:^  {e[P'(X,K^)]}2  (48) 

Thus  both  the  conventional  and  high  resolution  estimates  for  the 
frequency-wavenumber  spectrxim  are  consistent  and  can  be  made  to  be 
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unbiased.  Details  of  the  Complex  Wishart  distribution  and  derivations 
of  the  means  and  variances  for  the  estimates  are  given  In  [16],  [17], 

F.  Spherical  Wavefronts 

The  development  of  the  frequency-wavenumber  spectrum  has  been 
based  on  the  assumption  that  Incoming  signals  may  be  represented  by 
sums  of  plane  waves .  Since  plane  waves  are  represented  by  vector 
wavenumbers  the  frequency-wavenumber  spectrum  determines  wavefront 
direction  and,  for  a  given  temporal  frequency,  the  velocity  of  the 
wavefront.  In  many  applications  such  as  measurement  of  seismic 
events  from  great  distances  through  the  earth's  crust,  the  assump- 
tion of  plane  wavefronts  Is  a  good  one  Indeed.  For  electrophysiological 
data,  where  smaller  dimensions  are  found,  the  assumption  of  plane  wave- 
fronts  is  probably  less  valid;  however,  such  an  assumption  may  not 
invalidate  the  usefulness  of  the  computed  frequency-wavenumber  spectrum. 

An  analysis  of  the  frequency-wavenumber  method  applied  to  spher- 
ical wavefronts  indicated  that  velocity  and  wavenumber  magnitude  are 
determined  but  wavefront  direction  has  little  meaning.  This  comes  from 
the  fact  that  the  wavenumber  for  a  spherically  sjrmmetrlc  wavefront  is 
a  scalar  quantity.   Thus  if  the  spherical  wavefront  is  given  by 

2TTl(k  r-f  t) 
^  o   o 

S(r)  =^ (49) 


where:        k  =  scalar  wavenumber  (cycles/cm) 

f  =  temporal  frequency  (Hz) 

X  =  2iTf 
o     o 
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then  it  Is  shown  in  Appendix  C  that  the  frequency-wavenumber  technique 
gives  a  delta  function  about  k  =  k  and  X  =  X  . 

°  CO 


P(X,k)  =  6(X-X^)  6(k-k^)  (50) 


where:        P(X,k)  =  frequency-wavenumber  spectrum 

From  the  results  given,  it  is  clear  that  the  frequency-wavenumber 
technique  does  not  determine  wavefront  direction  for  the  spherical  wave- 
front.  For  the  case  of  spherical  wavefronts,  the  objective  should  be 
to  detejrmine  coordinates  for  wave  sources  relative  to  the  array  of 
sensors.  Accomplishment  of  this  objective  requires  the  use  of  methods 
other  than  the  frequency-wavenumber  method. 

To  apply  the  frequency-wavenumber  method  to  the  analysis  of 
electrophysiological  data  it  is  assumed  that  the  signals  may  be 
represented  by  approximate  plane  wavefronts  resulting  from  sources  at 
distances  large  compared  to  the  array  size  or  by  a  stunmation  of  plane 
wavefronts. 

G.  Spatial  Pre-filtering 

In  order  to  reduce  aliasing  effects  it  is  important  to  know  the 
signal  bandwidths  (spatial  and  temporal)  to  establish  maximum  sample 
intervals.   It  has  been  established  that  most  EEG  temporal  frequencies 
are  below  50  Hz  with  most  of  the  energy  at  2  Hz  or  below  for  averaged 
visual  evoked  responses  [18] .  No  limits  have  been  established  for 
spatial  frequencies;  therefore,  some  exploratory  flexibility  is 
maintained  for  initial  experiments  with  spatial  signals.  Spatial 
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pre-filtering  to  limit  wavenumber  content  and  reduce  errors  due  to 
aliasing  in  the  spatial  frequency  domain  was  considered  for  applica- 
tion to  the  experimental  data.  The  maximum  wavenumber  content 
allowable  without  spectral  overlap  was  deteinnined  by  electrode 
spacing  for  the  various  experimental  data  runs.  Two  types  of 
experimental  data  were  used:   1.  human  visual-evoked  response  data 
recorded  previously  by  the  University  of  Florida  Visual  Sciences 
Laboratory,  and  2.  data  recorded  from  the  neocortex  of  rats. 
Minimum  electrode  spacing  for  the  visual-evoked  response  (VER) 
data  was  two  centimeters  and  had  been  determined  previously  by 
crosstalk  considerations  for  scalp  measurements  [19].  Minimum 
spacing  for  the  rat  data  was  determined  by  available  precision  for 
array  construction  to  be  one  millimeter. 

Digital  spatial  filters  for  pre-filtering  the  experimental 
data  were  not  used  for  two  reasons:   1.  restrictions  on  array 
geometry  (uniform  sampling  rates)  imposed  by  available  multi- 
dimensional filtering  techniques  would  make  application  of  the  filter 
to  the  VER  data  impossible  and  2.  the  current  state-of-the-art  for 
implementation  of  digital  spatial  pre-filters  is  such  that  their 
characteristics  are  not  well  understood. 

Rather  than  introduce  ambiguity  and  possible  distortion  Into 
the  experimental  data,  an  alternate  method  for  checking  wavenumber 
content  was  used.  The  method  involves  a  comparison  of  the  wavenumber 
spectra  for  a  given  data  set  which  has  been  computed  for  two  different 
spatial  sampling  rates.   Such  a  comparison  was  made  for  data  from  a 
rat  in  which  three  by  three  matrices  of  two  and  one  millimeter  grids 
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respectively  were  used  for  recording  the  data.   This  method  allows 
detection  of  aliasing  errors  which  would  indicate  that  a  higher  spatial 
sampling  rate  was  required. 

It  is  expected  that  incorporation  of  spatial  pre-filtering  into 
the  analysis  technique  would  eventually  be  desirable.  For  further 
information  on  digital  filtering  and  spatial  filtering  the  reader  is 
referred  to  [20]  and  [21]. 

H.   Cost  and  Computation  Requirements 

Many  factors  were  involved  in  the  determination  of  computation 
time  and  cost  for  Individual  computer  runs;  however,  some  factors 
appeared  to  dominate  over  others.  Computation  time  and  therefore  cost 
was  proportional  to  the  number  of  electrodes,  the  total  number  of 
wavenumber  points  in  the  output  spectrum,  and  the  total  number  of 
san^le  points  of  the  spatio-temporal  field  (K  x  M  x  N).  For  the  band- 
pass computation  method,  increases  in  N  should  have  less  effect  because 
of  the  use  of  the  FFT  algorithm  in  computing  power  spectral  estimates 
of  the  N  temporal  data  samples.  This  was  verified  by  experimental 
results.  The  savings  in  computation  time  by  use  of  the  bandpass  method 
was  low  but  is  expected  to  increase  as  the  number  of  time  samples  is 
increased. 

Shown  in  Table  2  are  some  of  the  time  and  cost  figures  for 
various  experimental  runs.   These  values  could  be  reduced  by  the  use 
of  a  binary  program  deck  which  would  eliminate  compilation  time. 
This  was  not  done  for  the  experimental  investigation  because  of  frequent 
program  changes.  Also,  the  requirement  that  several  arrays  in  the  pro- 
gram must  be  dimensioned  exactly  make  the  use  of  an  object  deck  impractical 
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unless  a  fixed  electrode  geometry  were  to  be  used  in  analyzing  numerous 
experimental  data  trials. 

The  first  two  lines  in  Table  2  are  representative  of  the  two- 
dimensional  simulation  experiments.  The  VER  data  program  included  an 
interpolation  routine  which  added  to  the  required  computation  time. 


CHAPTER  V 
FREQUENCY-WAVENUMBER  SPECTRUM—SIMULATION 

The  validity  of  the  concepts  of  the  frequency-wavenumber 
analysis  technique  was  tested  via  a  series  of  simulation  experiments. 
Characteristics  of  the  frequency-wavenumber  spectral  estimate  were 
obp-erved  for  known  signal  and  noise  conditions.   The  effects  of 
certain  parametric  variations  were  demonstrated  by  these  experiments 
and  theoretical  concepts  were  verified  by  the  results. 

Because  of  its  relative  simplicity  and  lower  cost  to  implement, 
a  two-dimensional  analysis  was  used  for  those  experiments  in  which  the 
principle  to  be  demonstrated  did  not  depend  on  the  use  of  three- 
dimensional  computations.   Relative  cost  versus  performance  was  dis- 
cussed in  Chapter  IV.  The  simulation  experiments  consisted  of  software 
implementation  of  the  frequency-wavenumber  analysis  technique.  The 
computer  language  used  was  FORTRAN  IV  and  runs  were  made  on  an  IBM  360/65 
computer.  Listings  of  the  various  programs  used  are  given  in  Appendix  E. 
Storage  requirements  for  the  program  were  large  (120,000  locations  for 
a  two-dimensional  analysis  of  a  16-electrode  array  and  256,000  locations 
for  a  three-dimensional  analysis  of  a  l6^electrode  array);  therefore, 
implementation  on  a  small  computer  would  not  be  feasible  without  suitable 
program  segmentation.   For  fewer  electrodes  this  might  be  possible  and 
could  provide  useful  results.  However,  as  indicated  by  the  results  of 
simulation  experiments,  the  technique  worked  best  for  large  numbers  of 
electrodes . 
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A.  Array  and  Signal  Example  Parameters 

The  array  used  in  the  two-dimensional  simulation  experiments 
is  shown  In  Figure  9.   It  is  a  two-dimensional  rectangular  array  of 
sensors  with  equal  spacing  in  both  dimensions.   The  4x4  (16-electrode) 
array  was  used  for  most  of  the  experiments;  however,  the  2x2  (circled 
electrodes)  array  was  used  for  some  of  the  experiments  to  test  the 
resolving  power  of  the  array  for  limited  spatial  information. 

Shown  also  in  Figure  9  are  the  signal  parameters  used  in  the 
two-dimensional  simulation  experiments.  Note  that  the  spatial  and 
temporal  sampling  rates  were  conservatively  high  at  eight  times  the 
signal  frequency  and  wavenumber  instead  of  twice  the  maximum  frequency 
as  required  by  the  sampling  theorem.   For  analysis  of  practical  and 
non-artificial  data,  sampling  rates  slightly  higher  than  twice  the 
Nyqulst  rate  were  used  to  reduce  aliasing  effects  caused  by  non-ideal 
bandllmlted  and  wavenumber  limited  signals.  Application  of  the 
frequency-wavenumber  technique  to  electrophysiological  data  is  dis- 
cussed in  Chapter  VI.   It  is  emphasized  that  the  achievable  resolution 
for  a  given  number  of  sensors  and  a  given  number  of  data  points  is 
improved  for  sample  rates  nearer  to  twice  the  Nyqulst  rate  than  that 
achieved  with  a  sample  rate  of  eight  times  the  maximum  signal  frequencies. 

Computer  results  consisting  of  two-dimensional  plots  of  the 
frequency-wavenimiber  spectral  estimates  for  the  various  simulated  condi- 
tions are  Included  in  this  chapter  along  with  perspective  plots  and 
cross-sectional  views.   The  directions  shown  in  Figure  9  for  the  vector 
wavenumbers  were  selected  so  that  a  vector  drawn  from  the  peak  value  in 
the  frequency-wavenumber  spectral  plot  to  the  wavenumber  origin  represents 
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51  =  cos    [2TT(f^t  -   ^^x  -  Yj^y)] 

52  =  cos    [2TT(f  2t  -  32^  -  T2y)  3 

53  =  cos    [27r(f2t  -  33X  -  737)] 


COMPOSITE  SIGNALS 
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b)   SIGNAL  PARAMETERS 


Figure  9.   TWO-DIMENSIONAL  ARRAY  GEOMETRY  AND  PARAMETERS 
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the  wavefront  direction.   The  particular  convention  selected  is  not 
important  except  that  interpretation  of  results  must  be  consistent. 

Array  geometry  and  signal  parameters  used  in  the  three-dimensional 
simulation  experiment  are  discussed  in  Section  G.  Results  of  the  three- 
dimensional  simulation  are  given  in  the  form  of  a  series  of  two-dimensional 
wavenumber  plots  for  x  and  y  wavenumber  components  with  the  z  wavenumber 
component  being  constant  for  a  given  plot. 

B.   Comparison  of  High  Resolution  and  Conventional  Estimates 

For  the  first  experiment  a  comparison  was  made  of  the  obtainable 

resolution  by  the  high  resolution  estimation  procedure  relative  to  that 

obtainable  with  the  conventional  estimation  procedure. 

A  perspective  plot  on  a  two-dimensional  wavenumber  grid  is  shown 

in  Figure  10  for  results  obtained  by  the  two  methods.  The  peak  for  both 

methods  occurs  at  k  =  -.0625  and  k  =  +.0625  as  it  should.   It  is  clear 

X  y 

that  the  high  resolution  method  provides  a  sharper  peak  and  thus  improves 
the  resolution  over  the  conventional  method.  One  important  point  to  be 
emphasized  is  the  role  of  noise  in  the  high  resolution  estimate.  As  was 
discussed  in  Chapter  IV,  a  small  amount  of  incoherent  noise  is  added  to 
the  cross  power  spectral  density  matrix  to  insure  that  its  inverse  exists. 
The  level  of  the  injected  noise  has  an  effect  on  the  resolution  obtained 
in  the  high  resolution  estimate.  Thus  the  results  of  Figure  10  are 
noise-free  for  the  conventional  estimate  and  noisy  (10%  in  this  case) 
for  the  high  resolution  estimate.   It  is  shown  later  in  the  experimental 
results  that  the  high  resolution  method  can  be  improved  even  more  when 
the  noise  level  is  reduced. 
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CONVENTIONAL  FREQUENCY-WAVENUMBER  SPECTRAL  ESTIMATE 
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HIGH  RESOLUTION  FREQUENCY-WAVENUMBER  SPECTRAL  ESTIMATE 


Figure  10.      CONVENTIONAL  AND  HIGH  RESOLUTION  PERSPECTIVE  PLOTS 
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Shown  in  Figure  11  is  a  two-dimensional  plot  of  values  for  the 
conventional  frequency-wavenumber  spectral  estimate.   The  values  shown 
are  in  db  attenuation  so  that  the  peak  occurs  at  a  value  of  zero. 
Only  integer  values  are  printed  and  the  actual  spectral  values  are 
rounded  to  their  lowest  integer  values.   Figure  12  shows  a  plot  of  the 
high  resolution  estimate  under  identical  signal  conditions.   For  clarity 
values  of  attenuation  larger  than  20  db  are  not  printed. 

Figure  13  shows  a  section  view  of  the  results  obtained  by  the 
two  methods.  This  is  a  plot  of  spectral  values  along  the  k  =  .0625 
line.  "* 

C.  Array  Geometry 

From  the  sampling  theorem  it  can  be  shown  that  resolution  is 
inversely  proportional  to  the  data  length.  For  temporal  signals  the 
data  length  is  taken  to  be  the  record  length  in  seconds  and  the 
resultant  effect  is  on  frequency  resolution.   For  spatial  transforms 
the  data  length  is  equivalent  to  array  aperture  In  meters  and  the 
resultant  effect  is  on  wavenumber  resolution. 

An  experimental  run  was  made  to  test  the  effect  of  array 

aperture  on  wavenumber  resolution.   In  one  case  a  16-electrode  array 

was  used  with  an  aperture  of  3x  by  3y  (see  Figure  9).   In  the  other 

s  s 

case  a  4-electrode  array  with  an  aperture  of  x     by  y    ,   that  is,  one- 

s  s 

third  as   large  in  x  and  y  directions  as  the  16-electrode  array  was 
used.     A  section  view  of  the  results  obtained  from  the  high  resolution 
method  for  these  two  arrays  is  shown  in  Figure  14.     The  three-db  band- 
widths  are  indicated  and  it  is  observed  that   for  equal  noise  levels 
an  aperture  which  is  three  times  larger  yields  a  resolution  bandwidth 
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-j-  1.0   k  =  +.0625  cycles/cm 
16  electrodes 
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Figure  13.   CROSS-SECTION  OF  CONVENTIONAL  AND  HIGH 
RESOLUTION  SPECTRAL  ESTIMATES 
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Figure  14.   RELATIVE  EFFECT  OF  ARRAY  APERTURE  ON  RESOLUTION 
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which  is  one-third  as  large.  Thus  experimental  results  agree  with 
expected  results  as  implied  by  the  sampling  theorem. 

Additional  indication  of  the  resolution  dependence  on  aperture 
is  demonstrated  by  comparing  Figure  12  with  Figure  15  which  gives  the 
high  resolution  attenuation  plot  for  the  4-electrode  example. 

D.  Temporal  Window  Functions 

Two  temporal  window  functions  were  tested  to  determine  their 
relative  effects  on  the  frequency-wavenumber  spectral  estimate.  The 
rectangular  and  cosine  tapered  window  functions  used  in  the  simulation 
runs  are  shown  in  Figure  16.   The  temporal  data  segments  were  multiplied 
by  one  of  these  two  window  functions  before  Fourier  transformation  both 
for  the  discrete  frequency  computations  and  fast  Fourier  transform  (FFT) 
computations. 

For  the  discrete  frequency  computation  method  and  a  relatively 
high  Incoherent  noise  ratio  the  resolution  of  the  frequency-wavenumber 
spectrum  was  the  same  for  both  rectangular  and  cosine  tapered  temporal 
window  functions.  For  a  ratio  of  incoherent  noise  power  to  total 
signal  plus  noise  power  of  R  =  0.1  and  discrete  frequency  computation, 
the  frequency-wavenumber  spectrum  of  two  plane  wavefronts  was  identical 
for  both  the  rectangular  and  cosine  tapered  temporal  windows.   The 
spectrum  for  the  rectangular  window  is  shown  in  Figure  17  and  the 
spectrum  for  the  cosine  tapered  window  is  shown  in  Figure  18. 

As  the  incoherent  noise  ratio  is  reduced,  some  smearing  is 
introduced  into  the  two-wavefront  spectrum  for  both  rectangular  and 
cosine  tapered  temporal  windows.   Figure  19  shows  the  spectrum  for 
a  rectangular  window  with  an  incoherent  noise  ratio  of  R  =  10 
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RECTANGULAR  TEMPORAL  WI^^)OW 


1.0 


0       O.IT 


COSINE  TAPERED  TEMPORAL  WINDOW 


Figure  16.      EXAMPLE  TEMPORAL  WINDOW  FUNCTIONS 
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Figure  20  shows  the  spectrum  under  Identical  conditions  except  that 
a  cosine  tapered  window  function  was  used.  The  figures  show  that 
smearing  is  only  slightly  higher  for  the  cosine  tapered  window. 

The  above  results  apply  only  to  the  simulation  runs  where 
discrete  frequency  computations  were  used  to  compute  the  cross  power 
spectral  density.   For  the  bandpass  method,  the  results  are  different 
and  less  predictable.   It  is  shown  in  Section  F  that  for  the  bandpass 
spectral  computation,  the  cosine  tapered  window  function  causes  less 
error  due  to  sidelobe  interaction  and  produces  a  frequency-wavenumber 
spectrum  with  more  clearly  distinguishable  wavefront  directions. 
Further  quantitative  effects  of  incoherent  noise  ratio  on  resolution 
are  discussed  in  Section  E. 

E.  Additive  Incoherent  Noise 

The  relative  effect  of  additive  incoherent  noise  on  the 
resolution  obtained  with  the  high  resolution  frequency-wavenumber 
technique  was  tested  by  adding  known  amounts  of  noise  to  the  cross 
power  spectral  density  matrix  before  inversion.  For  the  example 
signals  used,  this  proved  to  be  a  necessary  step  to  insure  non- 
singularity  of  the  matrix.  Values  of  the  ratio  of  incoherent  noise 
power  to  total  signal  plus  noise  power  of  R  =  0.0001,  0.001,  0.01  and 
0.1  were  used  to  determine  the  effect  of  incoherent  noise  power  ratio 
on  resolution.  A  cross-sectional  view  of  spectra  obtained  for  a 
single  plane  wavefront  and  for  various  values  of  R  is  shown  in 
Figure  21.  ShoMi in  Figure  22  are  the  spectral  profiles  for  a  signal 
composed  of  two  plane  wavefronts. 
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Figure  21.   EFFECT  OF  INCOHERENT  NOISE  ON  SINGLE  WAVEFRONT 
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Figure  22.   EFFECT  OF  INCOHERENT  NOISE  ON  MULTIPLE  WAVEFRONTS 


6.6 


Ideally  it  is  desirable  to  have  zero  noise  level;  however,  the 
matrix  inversion  requires  in  many  cases  a  finite  noise  level  to  insure 
non-singularity  (  See  Appendix  F  ).  This  was  true  of  the  simple  plane 
wavefront  example  signals  used  in  the  simulation.  For  the  example 
signals  and  the  matrix  inversion  subroutine  used,  a  minimum  value  of 
R  =  6.0  X  10"   was  found  to  be  necessary  for  the  matrix  to  be  non- 
singular.  However,  because  of  the  spectral  smearing  at  low  noise  levels 
a  value  of  R  =  0.001  to  0.01  was  found  to  be  best.  For  data  recorded 
from  practical  sensor  arrays  it  is  expected  that  any  superficial  noise 
addition  would  not  be  necessary  because  of  the  inherent  noise  in  the 
data  and  its  general  complexity.  This  was  verified  by  results  obtained 
from  application  of  the  frequency-wavenumber  technique  to  the  analysis 
of  electrophysiological  data.  As  shown  in  Figures  21  and  22  resolution 
is  improved  for  both  single  and  multiple  component  signals  as  the  noise 
is  reduced.  However,  there  is  a  non-linear  interaction  present  which 
causes  the  resolution  to  improve  less  if  multiple  wavefronts  compose 
the  signal.  Figure  23  shows  the  high  resolution  db  attenuation  plot 
for  one  wavefront  incident  on  the  16-electrode  array.  The  incoherent 
noise  ratio  is  R  =  0.001.  The  improvement  in  resolution  for  reduced 
noise  level  can  be  seen  by  comparing  this  figure  with  Figure  12  where 
the  noise  ratio  was  R  =  0.1.  Figure  23  shows  that  even  the  closest 
wavenumber  points  to  the  actual  signal  values  had  spectral  magnitudes 
that  were  greater  than  20  db  down  from  the  peak  value.  Further  reductions 
in  noise  are  possible;  however,  noise  ratios  less  than  R  =  0.001  were 
found  to  introduce  smearing  into  the  spectral  estimate  in  the  form  of 
sidelobe  enhancement  when  more  than  one  wavefront  comprised  the  signal. 
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Spectral  plots  for  two-wavefront  signals  are  shown  in  Figures  24 
and  25  for  incoherent  noise  ratios  of  R  =  0.01  and  R  =  0.001  respectively. 
Cosine  tapered  window  functions  were  applied  to  the  temporal  data  before 
transformation.  As  was  mentioned  in  Section  D  the  spectral  smearing  seen 
in  Figures  24  and  25  was  only  slightly  greater  (within  1  db)  than  that 
observed  with  a  rectangular  temporal  window. 

F.   Bandpass  Signal  Spectral  Computation 

The  bandpass  spectral  computation  method  described  in  Chapter  IV 
was  simulated  by  applying  the  method  to  signals  representing  single 
plane  wavefronts. 

For  a  signal  represented  by  51  in  Figure  9  the  frequency-wavenumber 
spectrum  was  computed  using  the  bandpass  method  for  a  passband  from  six 
to  twenty  Hz.  The  spectrum  was  computed  for  both  a  rectangular  temporal 
window  and  a  cosine  tapered  temporal  window  with  respective  spectral 
plots  as  shown  in  Figures  26  and  27.   Spectral  sharpness  is  best  for  the 
rectangular  temporal  window  as  it  was  for  the  method  using  discrete 
computation  of  the  cross  power  spectral  density.   However,  for  signals 
containing  multiple  wavefronts  and  for  low  incoherent  noise  levels  the 
bandpass  spectrum  using  a  cosine  tapered  window  is  more  reliable  than 
the  bandpass  spectrum  using  a  rectangular  temporal  window.  This 
conclusion  is  based  on  simulation  experiments  for  the  two  window  func- 
tions in  which  the  two  wavefronts  were  more  distinguishable  for  the 
cosine  tapered  temporal  window.   This  might  be  expected  since  cosine 
tapering  is  a  form  of  smoothing  and  should  produce  a  spectrum  with  less 
error  due  to  sidelobe  effects.   An  unexplained  dependence  on  specific 
wavenumber  content  of  the  computed  spectrum  of  multiple-wavefront 
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signals  was  also  observed.   For  example  a  signal  composed  of  wavefronts 
SI  and  S3  (see  Figure  9)  was  successfully  resolved  by  the  frequency- 
wavenumber  method;  whereas,  a  signal  composed  of  wavefronts  SI  and  S2 
produced  an  ambiguous  spectrum  for  the  same  frequency  passband.   It  was 
found,  however,  that  judicious  selection  of  the  passband  effectively 
reduced  ambiguity  in  the  frequency-wavenumber  spectrum.  For  example, 
the  composite  signal  of  SI  (12.5  Hz)  and  S2  (15  Hz)  was  poorly  resolved 
by  the  frequency-wavenumber  spectrum  for  a  passband  of  six  to  twenty  Hz 
regardless  of  the  value  of  incoherent  noise  injected  into  the  signal  or 
the  choice  of  temporal  window.  The  same  example  was  run  with  a  different 
passband  (11.5  to  16  Hz)  which  produced  an  improved  estimate  of  the 
frequency-wavenumber  spectrum  over  that  obtained  with  the  larger  pass- 
band.  Figure  28  shows  the  narrow  bandpass  spectrum  for  a  rectangular 
window  and  Figure  29  shows  the  narrow  bandpass  spectrum  for  a  cosine 
tapered  window.  A  sharper  spectrum  is  achieved  with  the  rectangular 
window;  however,  a  more  accurate  spectrum  is  obtained  by  using  the  cosine 
tapered  window.  The  results  shown  in  Figures  28  and  29  are  consistent 
with  the  expected  result  of  cosine  tapering,  which  is  a  form  of  smoothing. 
In  general,  smoothing  produces  a  broader  spectral  estimate  with  higher 
statistical  reliability  and  less  error  due  to  sidelobes. 

Thus  the  results  of  experiments  using  the  bandpass  computation 
method  indicate  that  the  selected  ppssband  should  be  narrow  and  carefully 
selected  to  improve  the  frequency-wavenumber  spectral  estimate;  this 
may  call  for  spectrum  analysis  prior  to  data  processing.  Also  for 
improved  reliability  of  the  spectral  estimate,  a  cosine  tapered  window 
was  better  than  a  rectangular  window  when  using  the  bandpass  method. 
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Recall  from  Section  D  that  for  discrete  computation  of  the  cross  power 
spectrum,  both  rectangular  and  cosine  tapered  windows  produced  similar 
spectra, with  the  cosine  window  producing  slightly  (1  db)  more  smearing 
for  a  low  incoherent  noise  ratio. 

G.  Three-Dimensional  Simulation 

Because  of  the  complexity  and  projected  cost  of  three-dimensional 
frequency-wavenumber  computations  and  since  most  parametric  interactions 
were  sufficiently  demonstrated  by  a  two-dimensional  analysis  only  one 
simulation  experiment  was  designed  to  test  the  three-dimensional  technique. 

The  signal  consisted  of  a  single  three-dimensional  plane  wave- 
front  as  shown  in  Figure  30.  The  alght-electrode  cubic  array  shown  in 
Figure  30  was  used  in  the  analysis. 

Since  the  signal  was  a  single  plane  wavefront  and  by  reducing 
the  Incoherent  noise  ratio  to  R  =  10~  for  improved  resolution,  signif- 
icant savings  in  computation  time  and  cost  were  possible  by  selecting 
an  eleven-point  matrix  of  wavenumber  values  In  three  dimensions  instead 
of  the  twenty-one-point  matrices  used  in  the  two-dimensional  analyses. 

Results  of  the  three-dimensional  frequency-wavenumber  simulation 
experiment  are  shown  in  Figure  31.   Two-dimensional  attenuation  plots 
were  generated  for  an  x  and  y  wavenumber  grid  for  specific  values  of 
the  z-direction  wavenumber.  In  Figure  31a  is  shown  schematically  how 
the  spectral  plots  were  generated.   Figures  31b,  c  and  d  give  the 
attenuation  plots  for  z-dlrection  wavenumbers  corresponding  to  the 
signal  wavenumber  (center  plot)  and  one  wavenumber  point  each  side  of 
the  signal  wavenumber.   Because  of  the  high  resolution  obtained  for  the 
three-dimensional  simulation  all  spectral  values  were  greater  than 
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Electrode  spacing  =  2  cm. 


^y 


a)  ARRAY  GEOMETRY 


S  =  cos   [2ir(ft  -kx-ky-kz)] 

A  y  Z 


f  =  12.5  Hz 


k     =  -.075  cycle/ cm 

k  =  +.075  cycle/cm 

k  =  -.075  cycle/ cm 
z 


b)   SIGNAL  PARAMETERS 


Figure  30.   THREE-DIMENSIONAL  ARRAY  GECMETRY  AND  SIGNAL  PARAMETERS 
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twenty  db  down  from  the  peak  value  as  Indicated  in  Figure  31.  For  this 
reason,  only  the  three  attenuation  plots  shown  are  included. 

Thus  the  three-dimensional  frequency-wavenumber  spectral  computa- 
tion is  demonstrated.   For  the  application  of  this  technique  to  multiple- 
wavefront  signals  and  bandpass  computation  the  parametric  interactions 
discussed  earlier  should  have  analogous  effects.   Increases  in  array 
aperture,  number  of  data  segments,  and  number  of  data  points  which  are 
required  for  improved  resolution  and  statistical  reliability  produce  a 
greater  increase  in  computation  requirements  and  cost  for  the  three- 
dimensional  analysis  method. 


CHAPTER  VI 

INVESTIGATION  OF  ELECTROPHYSIOLOGICAL  DATA 
BY  FREQUENCY-WAVENUMBER  TECHNIQUE 

The  frequency-waveniunber  technique  was  applied  to  two  types 
of  electrophysiological  data  to  determine  if   the  technique  might 
provide  useful  information  about   the  data.      Human  visual-evoked 
response    (VER)  EEG  data  and  penicillin-induced  focal  epileptic 
activity  in  the  rat  neocortex  were  analyzed  via  the  frequency- 
wavenumber  method  to  demonstrate   the  application  of  the  technique 
to  these  data.      Succint  interpretation  of  the  results  appears 
premature  at  this  time;  however,    the  results   indicated  that  certain 
characteristics  of  the  data  might  be  categorized  by  application  of 
the  frequency-wavenumber  spectral  technique. 

A.     Electrophysiological  Signal  Characteristics 

For  the  data  analyzed,   a  power  spectral  analysis  was   first 
performed  to  determine  the  temporal  frequency  content.      Recall  from 
Chapter  V  that  spectral  smearing  and  ambiguity  in  the  wavenumber 
domain  were  less   for  those  examples  in  which  the  selected  temporal 
frequency  band  most   closely  approximated  the  true  frequency  content 
of  the  signal.     The  expected  wavenumber  content  of  EEG  signals  has 
not  yet  been  determined;  however,  physical  constraints  of  electrode 
size  and  array  spacing  determined  the  maximum  wavenumbers  which  pro- 
vided non-aliased  spectra.      Ideally  the  data  should  be  subjected  to 
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a  wavenumber-limited  spatial  pre-filtering  process  prior  to  application 
of  the  frequency-wavenumber  spectrum.  However,  a  more  expedient  method 
for  checking  the  wavenumber  content  was  designed  and  tested  using  the 
rat  data.   The  method  consists  of  computing  the  frequency-wavenumber 
spectrum  for  electrode  arrays  of  different  spatial  sampling  rates  and 
then  comparing  the  resultant  spectra  for  differences  due  to  aliasing. 

The  VER  data  had  been  previously  recorded  by  the  University 
of  Florida  Visual  Sciences  Laboratory  with  an  array  geometry  designed 
for  determining  characteristics  other  than  the  expected  wavenumber 
content.   For  a  description  of  the  experimental  method  for  recording 
the  VER  and  related  discussions  of  its  implications  see  [18],  [22] 
and  [23].   Section  B  gives  the  frequency-wavenumber  spectral  plot  for 
VER  data  obtained  from  a  subject  with  intermittent  exotropy  (subject 

DLH). 

For  data  recorded  from  the  rat  neocortex  the  spatial  sampling 
rate  was  restricted  by  the  techniques  available  for  the  electrode 
array  manufacture  and  by  the  electrode  diameter.  The  stainless  steel 
wire  electrodes  had  a  diameter  of  approximately  0.4  mm.  Thus  an  array 
of  these  electrodes  in  their  closest  spaced  arrangement  could  efficiently 
sample  signals  with  wavenumbers  no  higher  than  12.5  cycles /cm.   To  keep 
tolerance  errors  low  a  three-electrode  array  with  2,0  mm  spacing  was 
fabricated.   Sampling  intervals  of  1.0  mm  were  obtained  by  micro- 
positioning  of  the  array.   For  the  1.0  mm  sampling  interval  signals 
with  wavenumbers  as  high  as  5.0  cycles/cm  could  be  sampled  without 
aliasing  errors.   Spectra  produced  by  the  2.0  mm  and  1.0  mm  sampling 
grids  were  compared  to  determine  relative  errors  due  to  aliasing. 
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Section  C  describes  the  results  of  the  application  of  the  frequency- 
wavenumber  technique  to  the  analysis  of  focal  epileptic  discharges 
from  rat  neocortex.  Details  of  the  experimental  method  employed  for 
collection  and  analysis  of  these  data  are  presented  In  Appendix  D. 

B.  Analysis  of  Human  Visual-Evoked  Response 

Visual  evoked  response  data  recorded  from  an  array  of  sixteen 
electrodes  were  analyzed  via  the  frequency-wavenumber  method.   Individual 
signals  represented  averaged  EEG  responses  to  visual  stimulation  of  the 
subject  and  were  recorded  from  scalp  electrodes  positioned  over  the 
occipital  region  of  the  brain.  The  electrode  array  consisted  of  sixteen 
electrodes  located  on  two  concentric  circles  of  radii  2.0  cm  and  4.0  cm 
respectively.   Eight  electrodes  were  uniformly  spaced  on  each  circle. 
An  electrode  at  the  center  of  the  circles  was  used  as  a  reference. 
Electrical  response  activity  was  monitored  for  a  period  of  500  msec 
after  stimulus  presentation  and  the  average  of  fifty  such  responses 
was  used  as  the  VER  data  for  a  given  electrode  position.  All  VER 
responses  for  the  various  electrode  positions  were  referenced  to  the 
center  electrode  prior  to  analysis  using  the  frequency-wavenumber 
method. 

Analysis  of  the  VER  data  indicated  a  need  for  higher  temporal 
sampling  rates  than  required  by  the  sampling  theorem  for  the  expected 
maximum  temporal  frequency  content  of  the  signals.  This  conclusion 
was  based  on  the  requirement  that  each  time  record  be  divided  into  a 
number  of  segments  greater  than  or  equal  to  the  number  of  electrodes 
(16  for  the  VER  data)  and  on  the  observation  that  short  time  lapse 
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transient  events  comprised  the  visual-evoked  response.  Thus  a  large 
number  of  data  points  for  short  observation  times  required  an  increased 
sampling  rate.  Additional  temporal  sample  points  for  the  existing  VER 
data  were  obtained  by  a  linear  interpolation  between  the  true  sample 
values . 

An  initial  analysis  of  the  VER  data  using  the  frequency- 
wavenumber  method  was  performed  for  temporal  frequencies  of  2.0,  5.0 
and  10.0  Hz.  The  maximum  wavenumber  included  In  the  analysis  was 
given  by  the  inverse  of  four  times  the  electrode  spacing  or  0.125 
cycles/cm.  The  resultant  frequency-wavenumber  spectral  estimate 
for  the  500  msec  visual- evoked  response  is  shown  in  Figure  32.  The 
fact  that  all  spectral  values  were  at  either  zero  or  greater  than 
twenty  db  attenuation  levels  is  partially  explained  by  subsequent 
analyses  in  which  the  500  msec  response  time  is  divided  into  two 
250  msec  time  periods.  The  frequency-wavenumber  spectrum  for  the 
first  250  msec  is  shown  in  Figure  33  and  the  spectrum  for  the  last 
250  msec  is  shown  in  Figure  34.  These  computations  were  for  temporal 
frequencies  of  8.0,  9.0,  and  10.0  Hz  to  include  observed  power  spectral 
peaks  near  these  frequencies.   Inclusion  of  the  2.0  Hz  temporal  fre- 
quency In  the  spectral  computation  for  the  500  msec  record  length 
probably  contributed  to  the  ambiguity  observed  in  that  spectrum. 
Spectra  for  the  subdivided  response  indicate  that  spatial  signal 
sources  changed  with  time  for  the  500  msec  observation  period  so  that 
the  frequency-wavenumber  spectrum  for  the  total  observation  time 
exhibited  a  smearing  effect  of  different  transient  spatial  sources. 
This  observation  in  conjunction  with  the  required  large  number  of 
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data  points  for  input  to  the  frequency-wavenumber  spectral  estimation 
technique  formed  the  basis  for  the  conclusion  that  higher  temporal 
sampling  rates  are  required. 

As  indicated  by  Figure  33  the  frequency-wavenumber  spectrum  for 
the  first  250  msec  time  period  following  stimulus  presentation  gave  a 

peak  value  at  k  =  k  =0.   This  may  be  the  result  of  a  plane  wave 
X    y 

propagating  in  the  k  direction  or  may  be  the  result  of  a  synchronous 
active  discharge  exhibited  by  neuronal  populations  under  each  electrode 
and  triggered  by  a  common  source.  The  latter  explanation  is  considered 
most  likely  since  similar  results  were  obtained  with  signals  recorded 
from  the  rat  neocortex.   Details  of  the  result  of  application  of  the 
frequency-wavenumber  spectral  estimation  technique  to  analysis  of 
artificially  induced  focal  epilepsy  are  discussed  in  Section  C. 

The  frequency-wavenumber  spectrum  for  the  last  250  msec  of 
the  VER  exhibits  bands  of  uniform  maximum  and  minimum  values.  Again 
this  may  be  partially  explained  by  attributing  it  to  a  temporal 
smearing  of  numerous  disparate  spatial  sources.   If  the  problem  is 
analyzed  in  reverse  by  taking  the  inverse  spatial  Fourier  transform 
of  a  flat  spectral  band,  a  different  possible  explanation  is  obtained. 

Given  a  frequency-wavenumber  spectral  plot  as  shown  in  Figure 
35,  the  inverse  spatial  Fourier  transform  is  defined  by 


f(5f,X)  = 


=  1 


„/r  \\  -2TTik*x  ,5- 
P(k,X)  e       dk 

.00 

(51) 

W      W 

3  +— ^  y  +— ^ 
°  2    °  ^  -2TTi(k  X  +  k  v) 
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X    y   J  J 

■'        dydx 
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region  of  maximum  spectral 
values  (0  db) 


W 

I 


-^K 


Figure  35.   FREQUENCY-WAVENUMBER  SPECTRUM  WITH  MAXIMUM  VALUE  BANDS 


or 


-2Tri(3  X  +  Y  y)   sin  it  W  x  sin  tt  W  y 
f(x,A)    =e  °  °       ^ ^ 


TTX 


Tly 


(52) 


where:    P(k,A)  =  frequency-wavenumber  spectrum 
f(x,A)  =  cross  power  spectrum 

Inspection  of  (52)  indicates  that  a  frequency-wavenumber 
spectrum  with  flat  spectral  bands  may  be  represented  by  a  plane  wave 
modulated  by  a  two-dimensional  spatial  function.  Rewriting  (52)  as 


f(x,A)  =  A(x,y)  P(x,y) 


(53) 


where:    A(x,y)  = 


sin  IT  W  X  sin  tt  W  y 


P(x,y)  =  e 


TTX         Try 

-27ri(e  X  +  Y„y) 
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we  observe  that  A(x,y)  has  the  form  of  a  spherical  wave  modulating 
P(x,y)  which  has  the  form  of  a  plane  wave.  Alternately  if  equation 
(52)  iiS  written  as  shown  in  (54)  the  signal  may  be  represented  as 
two  separate  plane  waves,  each  of  which  is  modulated  by  a  one- 
dimensional  spatial  function. 


f(x,A)  = 


sin  TT  W  X     -2Tri3  x 
X  o 

e 


TTX 


sin  7T  W  y     -2-niy  y 


TTy 


(54) 


A  possible  explanation  is  postulated  for  the  flat  spectral  regions  in 

the  frequency-wavenumber  spectrum  of  the  VER  data  by  considering  the 

existence  of  active  neuronal  populations  within  the  recording  region 

which  modulate  wavefronts  propagating  across  the  region.  If  one 

accepts  this  hypothesis  then  wavefront  directions  as  given  by  3^  and 

Y  are  determined  by  the  centroid  of  flat  spectral  regions  in  the 
'o 

wavenumber  plane. 

In  any  case  it  is  demonstrated  that  the  frequency-wavenumber 
method  can  provide  additional  information  for  facilitating  the  under- 
standing of  visual  functioning  via  analysis  of  evoked  electrical 
activity. 

C.  Penicillin  Induced  Focal  Epilepsy  in  Rats 

The  exact  neuronal  mechanisms  underlying  the  spike-wave 
discharge  observed  in  the  EEG  of  epileptic  subjects  is  not  well 
understood  and  numerous  recent  research  efforts  have  been  directed 
toward  furthering  knowledge  of  the  causes  and  mechanisms  of  epilepsy. 
Focal  epileptogenesis  induced  by  freezing,  local  application  of 
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penicillin,  alumina  gel,  and  other  agents  has  been  studied  extensively 
in  frogs,  rats,  cats,  and  monkeys  [24],  [25],  [26],  [27]  and  [28].  The 
epileptogenic  focus  induced  with  local  application  of  penicillin  has 
been  shown  to  produce  spike  and  wave  activity  similar  to  human  post- 
traumatic epilepsy  [29]  and  [30].  This  fact  makes  the  study  of 
penicillin-induced  epilepsy  particularly  interesting  for  further 
investigation.  Associated  with  the  development  of  a  primary  penicillin- 
jnduced  epileptic  focus  is  the  development  of  a  secondary  or  "mirror 
focus."  Extensive  experiments  have  been  conducted  to  determine  the 
mechanism  and  pathways  by  which  a  mirror  focus  is  developed  [31],  [32] 
and  [33] .  Both  callosal  and  sub-cortical  pathways  have  been  postulated 
as  being  involved  in  the  development  and  maintenance  of  secondary 
epileptogenic  foci,  and  Isaacson  et^  al.  [33]  have  postulated  a  sub- 
cortical "pacemaker"  system  which  might  supply  synchronization  pulses 
which  trigger  both  primary  and  secondary  epileptogenic  discharges. 

Whatever  the  exact  mechanism  for  focal  transfer  it  seems  logical 
that  an  analysis  method  which  determines  wavefront  directions,  such  as 
the  frequency-wavenumber  spectral  estimation  method,  might  be  applied 
to  the  analysis  of  epileptiform  EEG  activity.  The  frequency-wavenumber 
method  was  applied  to  the  analysis  of  EEG  data  recorded  from  rats  with 
penicillin-induced  focal  epilepsy.  Details  of  the  experimental  method 
are  given  in  Appendix  D.  The  results  from  the  frequency-wavenumber 
analysis  of  data  from  three  animals  under  different  experimental  condi- 
tions are  included  in  this  section. 

Results  from  the  first  animal  were  significant  because  certain 
experimental  precautions  were  Indicated  as  being  important  specifically 
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for  application  of  the  frequency-wavenumber  analysis  technique  to  the 
experimental  data.  The  initial  experiment  used  a  four- elect rode  array 
of  twisted  bipolar  stainless  steel  electrodes  for  measuring  EEG  activity 
on  the  secondary  cortex.  Application  of  the  frequency-wavenumber  spectrum 
to  averaged  epileptogenic  discharges  from  the  electrodes  produced  a 
spectrum  which  indicated  the  presence  of  wavefronts  propagating  from 
directions  near  the  front  and  rear  of  the  rat  brain.  However,  subsequent 
observations  indicated  an  apparent  phase  reversal  for  the  average  wave- 
forms for  two  of  the  electrodes.  The  anfciguity  introduced  into  the 
spectrum  because  of  this  observation  emphasized  the  necessity  for 
maintaining  a  common  reference  for  all  electrodes.  Also,  for  the  four- 
electrode  array,  dimension  tolerances  were  poor  so  that  distortion  was 
introduced  into  the  frequency-wavenumber  spectrum  by  inaccurate  coordinate 
inputs.  Analysis  of  data  from  this  first  animal,  however,  did  confirm 
the  feasibility  of  applying  the  frequency-wavenumber  method  to  the  analysis 
of  averaged  epileptogenic  data.  Shown  in  Figure  36  is  the  frequency- 
wavenumber  spectrum  for  the  first  experimental  animal  for  temporal  fre- 
quencies of  2.7  and  4.7  Hz  and  with  phase  reversal  compensated. 
Referenced  to  the  coordinate  axes  in  Figure  36  the  front  of  the  rat  brain 
is  in  the  positive  k  direction.  The  array  was  located  on  the  right 
cortical  hemisphere  so  that  the  penicillin-injected  left  hemisphere  is 
in  the  negative  k  direction  for  the  frequency-wavenumber  plot.  As  can 
be  seen  from  Figure  36  the  wavefront  directions  were  no  longer  from  near 
front  and  rear  as  indicated  by  results  obtained  before  phase  reversal 
compensation.  Also  shown  in  Figure  36  is  the  indication  that  higher 
wavenumbers  may  constitute  the  signal.   Thus  the  need  for  a  higher 


93 


o      o 


.4      M      ^      M      >«      o      o 
o     o      o     o     o     o     o 


oooooooooo 

lijuiuiUiujaJUJUJiUUi  ujtuujujujujujiuiuuj 

ooaooooooa  oooocoosoo 

owomomoinoo  oouiaincmomo 


■a 

0) 

rt 
a 

f 

u 

0) 
01 

n 
p. 


N 

a 

4 
II 


cs 


O 

m 
en 


H 


(U 
M 

3 
tiO 
•H 


ooooeooooo 


I   I   I 


I   I 


94 


spatial  sampling  rate  Is  Implied.   For  data  from  animals  two  and  three 
a  common  reference  electrode  was  used  as  described  In  Appendix  D  and 
a  more  accurately  dimensioned  electrode  array  was  used  to  monitor 
cortical  electrical  activity. 

For  animal  number  two  the  electrode  array  was  positioned  to 
give  data  from  three  by  three  arrays  with  2  mm  and  1  mm  grids  respec- 
tively. A  power  spectral  analysis  of  the  data  indicated  a  spectral 
band  around  4  Hz.  Therefore  frequencies  of  3,  4,  and  5  Hz  were  used 
for  discrete  computation  in  the  frequency-wavenumber  spectral  technique. 
A  cosine  tapered  temporal  window  function  was  used.   Shown  in  Figure  37 
is  the  frequency-wavenumber  attenuation  plot  for  the  2  mm  grid  array, 
and  in  Figure  38  is  shown  the  same  spectrum  for  the  1  mm  grid  array. 
The  broader  spectral  peak  and  apparent  sidelobes  seen  in  Figure  37 
but  not  in  Figure  38  indicate  that  a  finer  sampling  grid  than  2  mm  is 
required  to  produce  an  alias-free  spectrum.  As  was  mentioned  for  the 
VER  data  analysis,  the  presence  of  a  spectral  peak  at  k^  =  k  =  0  may 
be  attributed  to  either  a  wavefront  propagating  in  the  z  direction  or 
a  synchronized  active  neuronal  discharge  across  the  array.   For  this 
animal  a  large  quantity  of  penicillin  powder  was  applied  to  the  primary 
cortex.which  resulted  in  very  strong  spike  discharges  on  both  primary 
and  secondary  cortical  surfaces.  The  high  level  spikes  completely 
masked  the  on-going  background  cortical  activity  so  that  some  informa- 
tion was  lost  in  the  data  used  for  analysis.   This  result  provided 
the  impetus  for  the  third  experiment  in  which  a  small,  carefully  applied, 
amount  of  penicillin  was  injected  into  the  rat's  left  cortical  hemisphere, 
Previous  experiments  have  indicated  that  the  quantity  of  penicillin 
inj6cted  is  one  factor  affecting  the  amplitude  and  rate  of  epileptic 
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discharge.  Also,  for  the  third  animal,  array  recordings  were  made 
on  the  primary  hemisphere  as  well  as  the  secondary  hemisphere. 

The  array  used  to  record  on  the  primary  hemisphere  had  a  2  mm 
grid ;  whereas ,  the  array  used  on  the  secondary  hemisphere  had  a  1  mm 
grid.-  A  power  spectral  analysis  indicated  strong  peaks  at  2  Hz  for 
both  the  primary  and  secondary  data  with  lower  peaks  at  7.8  and  13.7  Hz 
for  the  primary  hemisphere  and  lower  peaks- at  9.8  and  21.5  Hz  for  the 
secondary  hemisphere.   The  frequency-wavenumber  spectrum  for  a  temporal 
frequency  of  2  Hz  is  shown  in  Figure  39  for  the  primary  hemisphere  and 
in  Figure  40  for  the  secondary  hemisphere.   The  spectral  plots  for  the 
higher  temporal  frequencies  are  shown  in  Figure  41  for  the  primary 
hemisphere  and  in  Figure  42  for  the  secondary  hemisphere.   Inspection 
of  these  figures  indicates  very  small  differences  between  the  2  Hz  data 
and  the  higher  frequency  data.  The  spectra  for  the  mirror  hemisphere 
indicated  the  possibility  of  higher  wavenumber  content  and  the  need  for 
an  even  smaller  spatial  sampling  grid.  Increasing  spectral  values  in 
both  directions  along  the  k  axis  indicated  the  possible  existence  of 
spectral  peaks  in  those  directions  for  larger  wavenumbers.  The  most 
significant  difference  between  animals  two  and  three  is  that  the  lower 
level  spiking  produced  by  a  smaller  quantity  of  penicillin  used  for 
animal  three  produr.ed  a  frequency-wavenumber  spectrum  with  more 
complicated  features.   Thus  it  is  reasonable  to  state  that  lower 
spiking  levels  tend  to  obscure  less  of  the  other  EEG  activity  present 
in  the  averaged  waveforms.  Perhaps  the  spike  discharge  is  triggered 
by  a  much  weaker  signal  which  is  masked  by  the  spike  itself.   In  any 
case,  the  comparative  results  of  animals  two  and  three  have  indicated 
areas  for  further  research.  These  areas  are  discussed  in  Chapter  VII 
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which  summarizes  the  properties  of  the  frequency-wavenumber  spectrum 
as  applied  to  the  analysis  of  spatio-temporal  signals  and  specifically 
to  the  analysis  of  electrophysiological  data. 


/*V 


CHAPTER  VII 
DISCUSSION 

The  frequency-wavenumber  spectral  estimation  analysis  technique 
has  been  demonstrated  to  be  valid  theoretically  and  experimentally  when 
extended  to  include  analysis  of  three-dimensional  array  data  and  to 
include  data  for  finite  temporal  frequency  bands.  Parametric  analyses 
indicated  that  the  dominant  factors  influencing  wavenumber  resolution 
were  the  array  aperture  and  the  ratio  of  incoherent  noise  power  to 
total  signal  plus  noise  power.  Parameters  which  had  a  lesser  effect 
on  wavenumber  resolution  but  strongly  affected  computation  time,  ambi- 
guity, and  statistical  reliability  of  the  estimate  were  the  temporal 
window  function,  number  of  segments,  total  number  of  data  points  per 
sensor,  and  agreement  of  selected  temporal  frequencies  with  those 
comprising  the  signals.  Nominal  values  for  various  parameters  were 
established  which  provided  a  reasonably  economical  program  for  demon- 
strating the  properties  of  the  frequency-wavenumber  spectrum. 
Application  of  the  frequency-wavenumber  technique  to  the  analysis  of 
electrophysiological  data  was  demonstrated.  Results  of  these  analyses 
indicated  that  the  frequency-wavenumber  analysis  technique  may  provide 
information  useful  in  promoting  the  understanding  of  brain  function. 

Analysis  of  the  VER  data  indicated  the  need  for  a  temporal 
sampling  rate  higher  than  that  required  by  the  sampling  theorem  in 
order  to  resolve  short-time-lapse  transient  spatial  events  occurring 
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in  the  500  msec  time  period  immediately  following  presentation  of  the 
visual  stimulus.  The  concept  of  spatially  disparate  transient  sources 
contributing  to  the  observed  VER  has  been  discussed  by  Childers  [18] 
and  Mesa  [34].   Results  from  the  initial  frequency-wavenumber  analysis 
of  VER  data  tended  to  agree  with  this  concept.  The  concept  of  spatially 
disparate  sources  contributing  to  synchronous  discharges  is  further 
supported  by  Minsky  [35]. 

Frequency-wavenumber  analysis  of  penicillin-induced  epileptiform 
activity  recorded  from  the  rat  neocortex  indicated  that  spatial  sampling 
intervals  of  1  mm  or  less  are  required  to  reduce  aliasing  in  the  wave- 
number  domain.  The  spike  wave  discharge  observed  in  the  EEG  tended  to 
mask  other  cortical  activity.   Even  when  a  small  quantity  of  penicillin 
was  injected  in  the  cortex  to  produce  lower  level  spiking,  the  frequency- 
wavenumber  spectrum  still  indicated  a  dominant  effect  caused  by  the 
synchronous  spike  discharge.  Additional  experiments  in  which  the  spike 
was  removed  from  the  averaged  data  might  provide  more  information 
concerning  the  mechanism  by  which  the  spike  discharge  is  triggered. 
Previous  experiments  [33]  have  indicated  that  a  time  latency  of  10  to 
50  msec  exists  by  which  the  secondary  spike  lags  the  primary  spike. 
Results  of  the  frequency-wavenumber  analysis  indicated  that ,  although 
a  lag  may  exist  between  cortical  hemispheres,  each  hemisphere  tended 
to  exhibit  a  synchronous  discharge  across  the  cortical  surface  of  the 

hemisphere. 

Additional  application  of  the  frequency-wavenumber  analysis 
technique  to  electrophysiological  data  is  recommended  to  provide  a 
more  statistically  sound  basis  for  interpretation  of  the  results. 
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An  ejiperiment  using  the  closest  possible  array  spacing  should  be 
performed  to  determine  more  accurately  the  expected  wavenumber 
maximum  for  electrophysiological  signals.  Further  analysis  and 
implementation  of  wavenumber  limited  digital  filters  are 
recommended  to  insure  minimal  aliasing  error. 


APPENDICES 


APPENDIX  A 
CROSS  POWER  SPECTRAL  DENSITY 

Given  an  array  of  K  sensors  with  vector  positions  x . ;  j=l, ,K, 

a  covariance  matrix  is  defined  for  the  discrete  outputs  of  the  K  sensors 
as 

Pj^(mT,nT)   =  E[{s^(mT)   -   Sj}-{sJ(nT)  -  sj}]  (A-1) 

where:   E[»]  denotes  expected  value 
T  is  the  sample  time 
m,n  are  integers 

p  j,  is  the  j-ilth  element  of  a  K  x  K  covariance  matrix 
S.  and  Sg    are  mean  values 
*  denotes  complex  conjugate 

If  the  random  process  represented  by  the  sensor  outputs  is  considered 
to  be  wide  sense  stationary  with  zero  mean,  then  the  elements  of  the 
covariance  matrix  depend  only  on  the  time  difference  between  signals 
and  (A-1)  becomes 

p.^((m-n)T)  =  E[Sj(mT)  S*(nT)]  (A-2) 

The  cross  power  spectral  density  matrix  has  elements  defined  by 
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f.^(A)  =  F{p.^(m-n)}  =  I       p  j^(m)  e^"^  (A-3) 

where:     m  =  m  -  n  is  the  time  difference 

X  =  2TTfT  is  the  normalized  frequency 
f  =  frequency  in  Hz 
T  =  sample  period  in  seconds 
F{*}  denotes  Fourier  transformation 

From  Fourier  theory  the  coefficients  of  the  expansion  represented  by 

(A-3)  are  the  p.g  and  are  given  by  the  Fourier  coefficient  inversion: 
J*- 


Pjjl(m)  - 


£j,(X)  a-"  11  (A-4) 


This  method  for  conrputlng  the  cross  power  spectral  density  is  known 
as  the  correlation  method  and  consists  of  first  computing  the  covarlance 
(correlation)  matrix  and  then  taking  its  Fourier  transform.   Another 
method  known  as  the  direct  method  can  also  be  used  to  compute  the  cross 
power  spectrum.   The  major  differences  between  the  two  methods  lie  in 
their  implementation  as  related  to  computational  time  required  and  the 
statistical  reliability  of  the  resultant  estimate  of  the  cross  power 
spectrum.  The  direct  method  is  outlined  below. 

Given  a  discrete  signal  from  the  jth  sensor,  S.  CnT),  n  =  1, ,Nj 

we  can  define  its  Fourier  transform  as 


F.(X)  =  F{s.(nT)}  =  I  I      S  (nT)  e^'^^  (A-5) 

J        J         n=0  -■ 

An  alternate  definition  for  the  discrete  transform  which  is 
consistent  with  Fourier  theory  and  which  results  in  division  by  the 
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total  number  of  data  points  for  the  cross  power  spectral  estimate  is 
F^<^>  =472  I     S.(nT)  e^^^  (A-6) 

The  coefficients  obtained  by  the  use  of  (A-6)  differ  from  those 
obtained  by  the  use  of  (A-5)  only  by  a  scale  factor.  This  scale  factor 
must  be  accounted  for  in  the  definition  of  the  inverse  Fourier  transform. 
Since  the  purpose  of  the  frequency-wavenuraber  technique  is  to  detect 
relative  spectral  peaks  and  since  the  cross  power  spectrum  is  to  be 
normalized,  the  definition  (A-6)  is  used  for  discrete  Fourier  transform. 
Thomas  [36]  has  shown  that  a  consistent  estimate  for  the  power  spectrum 
is  obtained  from 


f  . (A)  =  lira  E  [F  (A)  F  (-X) ]  (A-7) 

for  a  single  sensor  output  S . .  This  may  be  extended  to  apply  to  the 
two-dimensional  cross  power  spectrum  by 


f.j(A)  =  lim  E  [F,(A)  F*(A)]      j,Jt=l, ,K     (A-8) 

Several  variations  on  the  direct  method  are  described  in  the  literature 
which  increase  the  statistical  reliability  of  the  power  spectral  estimate. 
One  such  technique  involves  segmenting  the  data  record  and  taking  Fourier 
transforms  of  each  segment.  The  spectral  estimate  is  then  averaged  over 
the  segments.  Both  overlapping  and  non-overlapping  segments  have  been 
used. 
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Capon  has  proposed  a  direct  segment  method  which  uses  non- 
overlapping  segments  for  estimating  the  cross  power  spectrum.  If  the 
Fourier  transform  of  the  mth  data  segment  from  the  jth  sensor  is  given 
by 

F,  (X)  =  470  I  a^  S.  (nT)  e^"^  (A-9) 

j  =  1, ,K 

ra  =  1, ,M 

N  =  number  of  data  points  per  segment 

Then  the  cross  power  spectral  density  function  is  a  K  x  K  matrix  whose 
elements  are  given  by 

fjt(«=S  ?/j„'«^L<"      J.W,— .K      (A-IO) 

•'  m=l 

The  cross  power  spectral  density  matrix  is  closely  related  to  the 
coherency  function;  however,  phase  information  is  lost  in  the  coherency 
function,  i.e. , 

2 
where:       Y.o(A)  =  coherency  function 

The  coherency  function  varies  from  zero  to  unity  and  is  a 
measure  of  the  linear  dependence  between  sensor  elements  [37]. 


APPENDIX  B 
SPECTRAL  REPRESENTATION  OF  SPATIO-TEMPORAL  FIELDS 

A  random  spatial  field  is  considered  to  be  homogeneous  if  the 
sensor  output  field  is  spatially  stationary.  That  is  the  elements  of 
the  cross  power  spectral  density  matrix  f^^CX)  for  fixed  X  is  a  func- 
tion of  the  vector  difference  x^  -  Xj^  between  sensors  j  and  Z  and  not 
the  vector  coordinates. 

Following  the  treatment  given  in  Yaglom  [12],  we  find  that  any 
homogeneous  random  field  has  a  spectral  representation  given  by 


Sj(nT)  = 


■^    f"  -i(nA+k'x,)  ^  .       . 

^  ^     Z(dX,dk)  <B-1) 


e 

-IT    -0° 


Vector 
where  k  is   the  vector  wavenumber  and  X  is  the  normalized  frequency. 
The  function  Z(AX,Ak)   is  a  random  interval   function  which  is  called 
the  "spectral  representation"  of  Sj(nT).      That  is,  the   function 
Z(A\,Ak)   associates  a  random  variable  with  each  interval   (AX.Ak)   and 
has  the  following  properties: 

1)  E[Z(AX,Ak)]  =  0  for  all  AA.Ak 

2)  Additivity 

Z(AX,A^k  +  A2k)  =  Z(AX,Aj^k)  +  Z(AX,A2k) 
if  A.k  and  A,k  are  disjoint  intervals  and 
Z(A^X  +  A^X.Ak)  =  Z(Aj^X,Ak)  +  ZCA^X.Ak) 
if  A^X  and  A2X  are  disjoint  intervals 
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3)  E[Z(A^X,A^k)  Z*(A2X,A2k)]  =  0 

if  A^A  and  A,^  are  disjoint  intervals  or 
if  A.ic  and  aS   are  disjoint  intervals 

4)  E[Iz(AA.Ak)|2]  =  P(A,k)  I^At 

Random  functions  satisfying  property  3)  are  called  random  functions 
with  uncorrelated  increments.  The  function  P(X,k)  is  the  signal 
power  density  contained  in  the  interval  (AX,Ak)  and  is  called  the 
frequency-wavenumber  power  spectral  density  function.  The  integral 
of  P(X,k)  over  all  interval  combinations  (AX.Ak)  gives  the  total 
signal  power  as  would  be  expected.  Thus  Z(AX,Ak)  may  be  considered 
as  analogous  to  the  voltage  spectrum  of  a  signal  and  P(X,k)  is  the 
power  spectral  representation. 

Using  (B-1)  and  the  definition  of  correlation  we  can  write  the 
cross  covariance  function  as 


p(m-n.x)  =  E[S.^  S^J 


(B~2) 


where:        p(m-n,x)  =  p^(m-n) 

for  a  homogeneous  field  (i.e.,  x  =  x.  -  x^)  or 


p(m-n,x)  =  E 


IT   00 


-i(mX+k*x.) 


Z(dX,dk) 


-TT  -0° 

Vector 


^  f"  +i(nX+k-x^) 


f( 


Z(dX,dk) 


-ir  -oo 
Vector 


(B-3) 
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Rearranging  gives 

^   r   -i[(m-n)X+k'(x.-  x„)]         .  -, 

(B-4) 


^  r   -i[(m-n)X+k-(x.-  x^)] 


p(m-n,x)  =     e  ^   ^   E[|  Z(dX,dk)  p] 


Vector 


or  on  replacing  the  time  difference  m  -  n  by  ra  and  recognizing  the 
definition  for  the  f requency-wavenumber  spectrum  we  have  finally 


IT  <» 


p(m,x)  =  j  I  e-it'"^-^^-^^  P(A.k)  f  dk  (B-5) 


-TT  -«> 

Vector 


If  both  sides  of  equation  (B-5)  are  multiplied  by  exp  (imX)  and  summed 
over  m,  we  have 


7T   00 


I       p(m.x)  e^"'^  =  f  f  e-^'^-^  P(X.k)  ^  dk  (B-6) 


IJ1=-00 


-IT  -" 

Vector 


or 


f(X,x)  = 


e-^^*^  P(X,k)  dk  (B-7) 

I 

Vector 


Recognizing  the  Fourier  transform  relationship  between  f(X,x)  and 
P(X,k)  in  the  spatial  domain  and  between  f(X,x)  and  p(m,x)  in  the  time 
domain,  we  can  write 
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P(X,k)  =   f  (X.x)  e^^*^  d3?  =  I 
J  m= 


00 


Vector  Vector 


Equation  (B-8)  shows  the  relationship  between  the  frequency-wavenumber 
spectrum  and  the  cross  power  spectral  density  function  or  the  cross 
covariance  function. 

To  show  how  the  frequency-wavenumber  spectrum  provides  informa- 
tion on  the  direction  of  an  incoming  wavefront,  consider  the  following 
example.   If  the  signal  consists  of  a  unity  magnitude,  monochromatic 
plane  wave,  then  the  signal  at  the  jth  sensor  is 

S.(mT)  =  exp  [-i(2Trf^mT  +  k^'x^)]  (B-9) 


where:        x.  is  the  sensor  vector  coordinate 
The  covariance  function  then  becomes 


p   (m-n,x.-  x^)  =  exp  [-i(2Trf^[m-n]T  +  k^'Cxj-  \)]  (B-10) 


or  if  ra  -  n  =  m  and  x.  -  x»  =  x  then. 


p.^(m,x)  =  exp  [-i(2Trf^mT  +  k^'x)]  (B-H) 


and  the  cross  power  spectral  density  function  becomes 

00 

f .^a.x)   =     I       exp   [-KZTTf^mT  +  k-x)  +  imX] 


in=-a> 


=  exp    L-ik  'x]      I       exp    [im(A-X  )] 


0  „ 


(B-12) 


where:  X^  =  2Trf^T 
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The  Identity 


I       exp    [27Tim(a-a  )]   =  6(a-a^)  (B-13) 


gives 


f(A,S)   =  2/T  exp   [-ik^'x]   6(A-X^)  (B-14) 


Substitution  of  this  expression  into  equation  (B-8)  gives  for  the 
frequency-wavenumber  spectrum 

oo 

P(X,k)  =  2tt  I  6(X-X^)  exp  [i(k-k^)-x]  dx  (B-15) 

_oo 

Vector 


or 


P(A,k)  =  6(A-X^)  6(k-k^)  (B-16) 


Thus  the  frequency-wavenumber  spectrum  allows  determination  of  the 
signal  frequency  and  wavenumber  and  specifies  the  plane  wave. 


APPENDIX  C 
FREQUENCY-WAVENUMBER  SPECTRUM  IN  SPHERICAL  COORDINATES 

Given  a  signal  in  spherical  coordinates,  S(r,6,(J)),  then  its 
Fourier  transform  is  given  by  [38] 

0°  TT  2Tr 
F(S,0,$)  =  J    J  S(r,6,<l)) 
0  0  0 

exp    [-2TTirS(cos  0  cos  6  +  sin  0  sin  6  cos   (<|)-$))] 
•   r^  sin  0  dr  d0  d(j)  (C-1) 

and 

00  71-  2ir 


S(r,e,(|))  =  I   I  \  ^(S,e,4>) 


0  0  0 

exp  [ZirirS  (cos  0  cos  6  +  sin  0  sin  8  cos  (({)-$))] 
•  S^  sin  0  dS  d0  d$  (C-2) 

If  S  is  a  function  of  the  spatial  coordinate  r  only,  then  a 

/' 

form  for  S  which  satisfies  the  wave  equation  is 

^^i^o'^   -2TTif  t 

S(t.r)  =  ^^-^—  e     °  (C-3) 

where:        k  =  scalar  wavenumber  (cycles /cm) 
o 

f  =  temporal  frequency  (Hz) 

X  =  2Trf 
o      o 
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Taking  the  temporal  transform  gives 


2Trik  r 
F(A,r)   =  ^ 6(A-y  (c-4) 


The  spatial  transform  for  r  dependence  only  is  then  given  by   [ 38] 

00 

P(A,S)  =  47r  J   F(A,r)   j^(2TrrS)   r^  dr  (c-5) 


where:  j^(2TrrS)  =  ^^2TrrS^'^   "  spherical  Bessel  function  of  order 

zero,  first  kind  (C-6) 

For  a  good  discussion  of  spherical  Bessel  function  see  [15], 

[39]  and  [40]. 

If  the  integration  in  (C-5)  is  completed  the  result  is  a  delta 

function  about  A  =  A  and  k  =  k  . 

o         o 


P(A,k)  =  6(A-A^)  6(k-k^)  (C-7) 


Thus  P(A,k)  is  the  frequency-wavenumber  spectrum  of  a  spherically 
sjnnmetric  wavefront  and  determines  the  frequency  and  scalar  wave- 
number. 


APPENDIX  D 
EXPERIMENTAL  METHOD 

The  method  used  for  obtaining  electrophysiological  data  from 
the  rat  neocortex  is  described.  The  experimental  animals  used  weie 
Long-Evans  hooded  rats  (male)  weighing  between  300  and  400  grams. 
After  being  anesthetized  with  i.p.  urethane  (ethyl  carbamate)  at 
a  level  of  11.6  ml  /kg  ,  the  animals  were  placed  in  a  stereo toxic 
frame  and  the  skull  was  opened  with  a  bone  rongeur  to  expose 
homologous  cortical  areas  on  both  hemispheres.  More  bone  was 
removed  from  the  hemisphere  which  was  to  be  monitored  with  an 
array  of  electrodes.  For  one  experiment  both  hemispheres  were 
exposed  for  array  monitoring. 

An  epileptogenic  focus  was  produced  on  the  left  cortical 
hemisphere  by  application  of  buffered  sodium  penicillin  G  powder. 
Application  was  made  with  the  aid  of  a  mounted  needle  with  a  blunt 
tip  which  contained  a  small  amount  of  penicillin  powder.  A  wire 
plunger  was  pressed  through  the  needle  to  a  level  which  barely 
punctured  the  dura  so  that  the  penicillin  powder  was  injected  sub- 
durally  into  the  primary  cortex.   For  some  animals  a  larger  amount 
of  penicillin  powder  was  placed  on  the  dura  to  produce  a  strong  and 
distinct  focal  epileptic  discharge. 

The  experimental  setup  for  monitoring  and  recording  the 
electrical  activity  of  the  brain  is  shown  in  Figure  D-1.  Enamel 
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Figure  D-1.  SCHEMATIC  DIAGRAM  FOR  MONITORING  AND 
RECORDING  ELECTRICAL  ACTIVITY  OF  RAT 
NEOCORTEX 
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insulated  stainless  steel  electrodes  with  bared  tips  were  used  to 
monitor  the  electrical  activity.  All  signals  were  bipolar  with 
respect  to  a  common  electrode  placed  on  the  primary  cortex.  Because 
of  equipment  limitations  only  four  electrode  positions  could  be 
recorded  at  any  given  time.  Therefore  a  precise  three-electrode 
array  was  designed  and  built  for  use  on  the  secondary  or  mirror 
cortex.  Thus  the  reference  electrode  and  three  array  electrodes 
were  monitored  and  recorded.  A  larger  array  was  simulated  by  micro- 
meter positioning  of  the  three-electrode  array  to  specific  locations. 
This  of  course  assumes  that  the  electrical  activity  is  statistically 
independent  of  the  time  origin;  that  is,  the  electrical  activity  of 
the  brain  may  be  considered  to  be  an  ergodic  random  process.  Because 
of  the  characteristic  spike  discharge  observed  in  the  electrical 
activity  of  a  brain  with  induced  focal  epilepsy  and  because  of  the 
inclusion  of  a  reference  electrode  signal  for  comparison  between  each 
successive  record  for  different  array  positions,  the  assumption  of 
time  origin  independence  is  considered  to  be  justified. 

Signals  from  the  electrodes  were  fed  directly  into  a  Grass 
Model  7  Polygraph.  A  low  frequency  cutoff  of  1  Hz  was  used.  An 
upper  cutoff  frequency  of  35  Hz  was  found  to  be  best  for  reducing 
noise  while  maintaining  good  reproduction  of  the  spike  discharges. 
The  J6  outputs  of  the  polygraph  were  calibrated  using  the  internal 
calibration  source  of  the  model  7  to  give  a  1  volt  rms  (2.8  v.  p-p) 
output  as  measured  on  an  oscilloscope.  The  J6  outputs  were  fed 
directly  into  an  FM  multiplexer  and  then  recorded  on  one  channel  of 
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a  stereo  tape  recorder.  The  other  stereo  channel  was  used  for  voice 
recording.  Records  of  approximately  five  minutes'  length  were  obtained 
for  each  electrode  position. 

Signals  from  the  various  electrodes  were  later  fed  into  a 
Digital  Equipment  Corporation  Lab-8  signal  averaging  system.   (See 
Figure  D-2.)  This  system  includes  a  precision  A/D  converter  which 
in  conjunction  with  the  Advanced  Averager  software  provides  a  means 
for  averaging  a  specified  number  of  epileptic  spikes.  Periods  of 
time  before  and  after  the  spike  were  included  in  the  displayed 
average  by  specifying  either  negative  or  positive  trigger  delays. 
Synchronization  pulses  were  produced  by  a  level-sensitive  Schmltt 
trigger  in  the  Lab-8.  The  averager  is  triggered  by  an  incoming 
signal  which  first  exceeds  the  upper  threshold  of  0  V  and  then  exceeds 
in  a  negative  direction  the  lower  threshold  (adjustable  from  -0.5  to 
-2.5  V).   The  threshold  adjustment  in  conjunction  with  the  gain 
adjustment  of  the  trigger  amplifier  was  used  to  discriminate  between 
spikes  and  background  activity.  All  averages  were  triggered  on  the 
signal  recorded  between  the  primary  and  reference  electrodes  which 
remained  in  the  same  positions  throughout  the  experiment.  This 
signal  was  denoted  the  reference  signal. 

Parameters  of  the  Advanced  Averager  were  specified  to  provide 
averages  of  50  sweeps  of  256  data  points  sampled  at  2  msec  to  4  msec 
sample  times  for  sweep  durations  of  512  msec  and  1.024  second  respec- 
tively.  Completed  average  waveforms  were  recorded  on  Polaroid  film 
and  digital  paper  tape  records  were  punched  for  each  waveform. 
These  tapes  were  used  as  inputs  to  an  8K  Fortran  program  which 


122 


FM  MULTIPLEXER 


AMP 


LAB-8  ADVANCED 
AVERAGER 


1 


SCHMITT 
TRIGGER 


2  ANALOG  INPUTS 


DIGITAL  TAPE 


PRINT 


DC  REMOVAL 


KEYPUNCH 


DB 

ATTENUATION 
PLOT 


IBM  360 


FREQUENCY- 
WAVENUMBER 
SPECTRAL 
ANALYSIS 


POWER 
SPECTRAL 
ANALYSIS 

— ir~ 


PUNCH 


J 


Figure  D-2.  PREPARATION  OF  RAT  DATA  FOR  FREQUENCY-WAVENUMBER  ANALYSIS 
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removed  the  DC  component  and  printed  the  data  values  In  a  suitable 
form  for  keypunching.  Because  of  amplitude  variations  in  the  signal 
from  the  reference  electrode  over  a  period  of  time  and  because  only 
three  electrode  signals  in  addition  to  the  reference  signal  were 
available  for  averaging  at  any  given  time  it  was  necessary  to  ensure 
that  all  reference  spikes  had  the  same  time  origin.  Visual  inspection 
of  printed  digital  values  for  the  reference  signals  made  it  possible  to 
specify  delay  units  required  for  ensuring  time  coincidence  of  the 
reference  signals.  The  delay  units  were  then  added  to  all  printouts 
for  averages  obtained  with  a  given  reference  signal  before  printout. 

A  power  spectral  analysis  was  performed  using  digital  tapes 
with  DC  components  removed  as  inputs  to  an  8K  Fortran  program.  The 
spectral  analysis  program  utilized  a  ten-percent  cosine  tapered  window 
function  and  an  FFT  algorithm  [41]  to  determine  spectral  content  of 
the  averaged  waveforms.  Using  data  cards  keypunched  from  the  data 
printout  and  selecting  temporal  frequencies  based  on  the  power 
spectral  analysis  the  frequency-wavenumber  technique  was  employed  to 
determine  wavenumber  content  of  the  focal  epileptic  discharges. 


APPENDIX  E 
PROGRAM  LISTINGS 
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APPENDIX  F 
SINGULARITY  OF  THE  CROSS  POWER  SPECTRAL  DENSITY  MATRIX 

For  reasonable  data  comprising  the  signal,  the  cross  power 
spectral  density  matrix  has  an  inverse  if  the  condition  M  ^  K  Is 
satisfied.  However,  under  certain  conditions,  the  signal  spatial 
characteristics  are  such  that  the  matrix  is  singular  unless  a  small 
amount  of  incoherent  noise  is  added  to  the  matrix  before  attempting 
to  invert  it.  A  cross  power  spectral  density  matrix  can  be  shown 
to  be  singular  if  there  is  a  linear  dependence  among  any  two  of  its 
rows  or  columns.  A  simple  example  demonstrating  linear  dependence 
is  the  matrix  obtained  for  a  plane  wavefront  at  normal  incidence 
to  a  linear  array.  A  more  complicated  example  is  given  by  the 
following  analysis. 

Consider  the  data  for  the  jth  sensor,  mth  segment  for  a  single 
plane  wavefront  signal: 

S.^(n)  =  exp  [27ri{(n-l  +  N(m-l))foTg  -  k^'x J]         (F-1) 

Using  the  block  averaging  method  for  computing  the  cross  power 
spectral  density  matrix,  we  have 
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The  elements  of  the  cross  power  spectral  density  matrix  given  by 
(F-2)  for  a  plane  wavefront  are  pure  exponentials.  Since  the  matrix 
is  Hermitian  and  has  pure  exponential  elements  it  can  be  shown 
that  the  matrix  is  singular.  This  is  illustrated  by  computing  the 
determinant  for  a  simple  2x2  matrix. 
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Thus  an  example  is  given  of  a  singular,  non-negative  definite 
matrix  for  a  deterministic  signal.  Singularity  of  the  cross  power 
spectral  density  matrix  for  more  complicated  deterministic  signals 
is  not  proven  in  general  and  may  be  dependent  upon  ma;:;nitudes  of  the 
wavenumber  and  of  the  electrode  spacing. 

A  method  for  making  the  cross  power  spectral  density  matrix 
positive  definite  and  thus  insuring  that  its  Inverse  exists  is 
described  In  Section  D  of  Chapter  IV. 
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